Fractional models in electrical impedance spectroscopy data for glucose detection

Biomedical Signal Processing and Control 40 (2018) 180–191

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Fractional models in electrical impedance spectroscopy data for glucose detection夽 Oscar Olarte a,∗ , Kurt Barbé b a Department of Mechanical Engineering, Acoustic and Vibration Research Group (AVRG) at Vrije Universiteit Brussel (VUB), Pleinlaan 2, 1050 Elsene, Belgium b Department Mathematics (DWIS), Faculty of Sciences, Vrije Universiteit Brussel (VUB), Pleinlaan 2, 1050 Elsene, Belgium

a r t i c l e

i n f o

Article history: Received 22 November 2016 Received in revised form 25 August 2017 Accepted 16 September 2017 Keywords: Fractional models Modeling Odd random phase excitation Sensor Glucose Electrical impedance spectroscopy

a b s t r a c t The article presents a methodology to discriminate glucose levels using electrical impedance spectroscopy technology. The method is based on an adequate estimation and assessment of the impedance data followed by identification of a general rational fractional model. The methodology is illustrated on a group of saline–protein–glucose solutions at physiological concentrations, and shows the ability of the fractional models to discriminate glucose levels. The method exhibit significant differences in the zero position of the fractional model for different glucose concentrations allowing discriminate the glucose effect on the impedance data along different matrix solutions. The results based on fractional model method are compared with classical circuit models and rational integer models. Even when the different methodologies perceive variations in the impedance data given glucose alterations, the fractional models present low uncertainty allowing discriminating small glucose alterations in the physiological range. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Conditions for a continuous blood glucose monitoring include high accuracy, high sensitivity and short measurement time [1–4]. These conditions establish electrical impedance spectroscopy (EIS) as one of the most suited techniques, since it allows distinguishing, in a single measurement, multiple individual contributions from different constituents of the system under test [5–10]. Unfortunately, the relation between the structure of the system, the analyte concentration and the impedance behavior is not straightforward. The relation can be complex and requires a suitable model that allows interpreting the measurement results while reduces the uncertainty in parameters. There are different approaches to establish a system’s model. The ideal method is to develop a mathematical model that correlates and quantifies each one of the physical–chemical processes that probably occur in the system. However, such a method requires plenty knowledge and critical understanding of all processes that could take place in the system. Other extended alternative is to use analog electrical circuit models where the circuit’s elements and connections represent

夽 This document is a collaborative effort. This work is financially supported by SRP-OPTIMech Vrije Universiteit Brussels (VUB). ∗ Corresponding author. E-mail addresses: [email protected] (O. Olarte), [email protected] (K. Barbé). http://dx.doi.org/10.1016/j.bspc.2017.09.017 1746-8094/© 2017 Elsevier Ltd. All rights reserved.

physical processes that are likely present in the system. In this sense, the designed circuit is believed to be the most appropriate design [5,6]. However, there exist difficulties to obtain the exact or correct circuit model for a particular process or system. The principal ones can be reported as: (1) it is possible to find more than one equivalent circuit that follows the spectrum impedance. (2) From the different circuits that can be identified, it is not possible to say what configuration is correct, and (3) it is always possible to obtain a more accurate fitting over the experimental data adding more circuits elements [11]. Fractional models for NaCl-glucose solutions √ have been previously reported in [12] under the fix s or Warburg variable showing important result reducing the model order. In the current work a general fractional model is introduced and assessed over different glucose solutions showing the capability of the method to discriminate glucose content at physiological levels. In order to develop a glucose measurement system based on EIS technology, the present article is concerned in an appropriate excitation signal strategy, followed by data analysis and a convenient modeling procedure. In general two steps are considered: • Measure and quantify the quality of the measured data: A convenient excitation signal that improves the EIS experimentation is a multisine excitation. This signal reduces the measurement time, offers full control of the excitation harmonics and shape of the power spectra keeping the properties of pseudorandom signals. This excitation allows estimating the additive noise and discriminate between even and odd nonlinear contributions present in

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

the system, while allowing to calculate the best linear approximation (BLA) of the impedance, as well as assessing the stationarity of the process [13]. All this information is important since electrochemical systems are prone to be highly non-linear, while out of controlled conditions could present high levels of noise. In this sense, the quality of the measurements is quantified based on the level of noise and levels of non-linear distortions present in the system. • The second step is identification based on fractional models. Fractional order systems are described by integro-differential equations involving fractional order derivatives. They are a generalization of the traditional calculus. Similar concepts and tools are applied, but with a wider generalization. [14–17]. The drawback, however, is that fractional models are non-linear in the parameters and require an initial estimation to guarantee convergence to the solution. In the present method the initial estimation is obtained by an algorithm based on the linear integer model in the Laplace domain [18]. The methodology is illustrated on a group of fundamental solutions (saline–glucose and albumin–saline–glucose) where the electrolyte concentration have been modified from average physiological to hiponatremia and hipocloremia levels (low levels of Cl and Na). In all cases the discrimination of the glucose was possible using fractional models, while circuit models and integer models showed high uncertainty avoiding discriminating the glucose levels. In this sense, the proposed methodology presents high capabilities to improve analyte detection based on EIS. There are several areas in science and engineering where fractional models have been successfully applied. Some examples in medical/biomedical applications are: fractional model structure describing the hemodynamics in functional magnetic resonance imaging (fMRI) [19], fractional order models to quantify the respiratory mechanical properties of the respiratory system [20], description of the fractional behavior of the deoxyribonucleic acid code [21], fractional models to understand the dynamics of human immunodeficiency virus infecting the CD4+T lymphocytes [22,23], fractional calculus models of complex dynamics in biological tissues [24]. In the field of electrochemistry, the use of fractional elements to model the impedance spectral data are described in [6,25]. In [26] the data analysis of mass-transfer electrochemical systems using fractional models is analyzed showing reduction in the number of parameters compared with integer models. A survey using fractional circuits models to fit impedance data in biomedicine and biology are presented in [27]. A review of the history of fractional calculus is found in [28]. The article is organized as follows: Section 2 the methodology for glucose analyte detection based on EIS is presented. In this section the excitation signal and the fractional model identification procedure are described. Section 3 the result and analysis of the impedance measurements and modeling are presented. Section 4 compares the results using fractional models with those obtained by circuits models and rational integer models. Section 5 presents the general conclusions of the present work. 2. Methodology The methodology includes two general steps: Estimates the BLA of the impedance system, assess the measured data, and identify a parametric fractional model. 2.1. Impedance measurement From the different signals that allow broadband excitation (maximum length binary sequence, noise excitation, chirp, discrete

181

interval binary sequence, multisine excitations among others), the random phase multisine reduces the measurement time, allows full control of the excitation spectra while allows identifying the level of noise and the level of non-linearities present in the system. The random phase multisine signal is a periodic signal formed by the sum of N harmonically related sine waves with amplitudes Ak and uniform distributed random phases k in the interval [0, 2[ 1  r(t) = √ Ak sin(jωk t + k ) N N

(1)

k=1

ωk are the excited frequencies (ωk = 2k f ) and fs the sampling N s frequency. From the different signals that could be generated by (1) this article uses an odd random phase multisine (ORPM) signal. In this signal most of the odd harmonics are excited, called excited odd harmonics. The even harmonics are not exited, called non-excited even harmonics. And in a group of predefined and consecutive odd harmonics one of them is randomly non-excited (Ak = 0), called nonexcited odd harmonics. In the multisine the lowest frequency is the fundamental and define the period of the signal. The maximum frequency is limited by the Nyquist theorem. This kind of signal is periodic, and from here the leakage is not a problem in the impedance estimation. In order to have good impedance estimation the measurements need to be done in steady state. To reach this condition a waiting or stabilization time (tW ) is required. This time depends on the time constant of the system and the frequency resolution [29,13]. For P measured periods, the measurement time using ORPM excitation (tORPM ) is: tORPM = tW + P

1 fmin

(2)

with (fmin ) the lower excited frequency. The waiting time (tW ) is considered only once, which reduces the measurement time compared with (equivalent) single-sine excitations, where the waiting time (tW ) need to be counted for each excitation frequency. On the other hand the ORPM is a broad band excitation signal and reduces the measurement time compared to swept–sine excitations [13,30]. The reduction of the excitation time is important since diminish de risk of estimating the impedance during changes in the system. However, it does not imply that the system behaves stationary. In order to assess the presence of non-stationary contributions, the levels of noise at excited and non-excited frequencies can be used. The variance of the noise should be a smooth function of the frequency. So, the noise level at excited and non-excited frequencies should be approximate the same. If the system is non-stationary, the noise estimation at the excited frequencies will contain information not only of the noise but of the non-stationary behavior in the system, generating a difference between the noise estimations (at excited and non-excited lines) [31–34]. More information on the measurement analysis using multisine excitations can be found in [35,36]. 2.2. Identification based on fractional models Once the measurements have been assessed and the BLA of the impedance system estimated, the next step is to establish the relationship between changes in the impedance given changes in the system by variations in glucose concentration. In order to establish such relation fractional models are developed. Modeling a fractional system using integer order models require a high number of parameters to reach a satisfactory approximation [37]. The high number of parameters will increase their uncertainty making difficult to identify small changes in the impedance

182

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

or |z1 | ≤ |p1 | ≤ |z2 | ≤ |p2 | ≤ · · · ≤ |pm−1 | ≤ |zm |

(8)

Each one of these group of roots can be reduced to a fractional pole or zero. The poles or zeros that do not follow the inequalities (7) or (8) will remain as an integer pole or zero (GR (ω) in (5)). An approximation of the fractional system that crosses through the staircase generated by the integer poles and zeros can be described by the compression formulas (9) to (11), where ˛comp represent the power of the fractional pole/zero, ωcomp the location of the fractional pole/zero and Kcomp the associated gain. The values ωc (k) are the locations of the roots in (7) or (8). Details of the proof and the fractional compression concept can be found in [18]. Similar approaches can be found in [37,38]: Fig. 1. Concept of the fractional order system approximation by an integer system.

as expected from glucose at physiological levels. In order to illustrate the high number of parameters needed by an integer model to simulate a fractional system, we introduce the single fractional pole G(s) =

1

(3)

(s − pi )˛

This real-value fractional pole presents a slope of −20˛ dB/decade in a Bode diagram, while an integer pole generates strict slopes of −20 dB/decade. In this sense, if we want to approximate the fractional system by an integer one, the approximation will be done by a group of integer poles and zeros generating alternate slopes of ±20 dB/decade, as is illustrated in Fig. 1. From here it is clear that the best approximation of the fractional system by an integer one, will requires a large number of intercalated integer poles and zeros: G(s) =

M−1

1 (s − pi )

˛

i=0 M M→∞

= lim

i=0

(s − zi )

(4)

(s − pi )

Note that the number of poles and zeros can be high as a group of poles and zeros are required to obtain a proper fit of the data. The increase model complexity implies a higher estimation variance while a fractional model improves the model variance as well as potential interpretation of its parameters. Fractional mathematics reduce the complexity of the integer models, which is important since the information should be distributed along the parameters. The reduction of the number of parameters and model complexity allows a direct understanding of the model and the system that produces such behavior. In general, a system can be modeled by



Gfrac (ω) = K



M 

m=1

1 (jω − pfrac (m))

˛m

GR (ω)

(5)

where M is the number of fractional poles/zeros, ˛m is the fractional multiplicities defined between 0 and 1. The multiplicities (˛m ) can be negative which will correspond to a fractional zero. The function GR (ω) refers to the integer part of the system and K the gain. In order to have a first estimation of the location of the fractional system, the method starts from the integer model system

mB

m=1 G(s) = K m A

m=1

(s − zm ) (s − pm )

(6)

The procedure of reducing the integer system (6) to a fractional one starts identifying the groups of poles and zeros that present intercalated locations as |p1 | ≤ |z1 | ≤ |p2 | ≤ |z2 | ≤ · · · ≤ |zm−1 | ≤ |pm |

(7)





1  log (ωc (2k + 1))/(ωc (2k + 2))   m−1 log (ωc (2k + 1))/(ωc (2k + 3)) m−2

˛comp =

(9)

k=0



1  = ωc (2k + 1) m−1 m−2

ωcomp

k=0

ωc (2k + 2) ωc (2k + 3)

G(0)  ˛+1 ωc (2k + 1)ωc˛−1 (2k + 2) m−1

(10)

m−2

Kcomp =

(11)

k=0

Once the fractional model structure and initial values have been calculated, an non-linear optimization process is carried out minimizing the cost function (12) with initial values given by the compression procedure equations:



N/2−1

arg min ˛m ,pfrac(m),K



2

Wk ∗ GBLA (ωk ) − Gfrac (ωk )

(12)

k=1

where GBLA (ω) in the best linear approximation of the impedance data, Gfrac (ω) is the fractional model defined in (5) and the weight factor, Wk , is the inverse of the variance of the impedance measurements calculated along P periods of the multisine excitation. 2.3. Considerations The proposed methodology considers two general steps: an appropriate measurement and quantification of the quality of the recorded data, followed by a convenient modeling procedure. The identified model contains a reduced number of parameters which reduces the uncertainty which allows discriminating small impedance changes. The proposed methodology is general and can be applied with no restrictions under the framework of linear time invariant systems. Any process out of this assumptions will require practical and experimental actions, or specific models for the time variation and nonlinear behavior. The design of the ORPM signal requires the selection of the minimum and maximum excitation frequencies, the resolution and the amplitude. This selection needs to be done knowing the system since, it could alter the linearity and the time invariance assumptions. A high amplitude excitation could generate or activate nonlinearities, while a long acquisition time could result in a time variant system (if the system changes faster than the recording procedure). Nonlinear behavior can be modeled by nonlinear models. These models are in general complex and generate high number of parameters. Nonlinear models are out of the scope of the current work and the reader can find information about this topic in [39–41]. A practical and well known alternative to keep the model in a linear state, or reducing the nonlinear effects, is to reduce the excitation amplitude. However, the price to pay is low SNR.

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

183

Table 1 Methodology procedure. Excitation Estimation Assessment

Modeling

Design the ORPM signal as in (1). Select the frequency resolution f0 , the frequency band, the number of periods, the excited and non-excited harmonics. Estimate the Impedance, noise and non-linear distortion levels. A general description can be found in Section 2.1 and a extended analysis in [13,29,34–36] * Asses Linearity → if the system is linear: the noise level and the total standard deviation are at the same level. → if the system is non-linear: reduce amplitude excitation to reduce non-linear behavor. * Asses Time in-variance. → if the system is time in-variant: soft transition between the noise level in the excited and non-excited harmonics is present. → if the system not time in-variant: reduce the measurement time (reducing the number of periods and/or increasing the lower excited frequency). Define fractional model as in (5). Identify the integer model as in (6). Identify the poles and zeros chains as in (7) and/or (8). Calculate the initial values of the fractional model (5) by Eqs. (9), (10) and (11). Define final fractional model by minimizing the cost function (12).

Table 2 Solutions.

Solution 1 Solution 2 Solution 3

NaCl [mg/dL]

Albumin [g/dL]

Glucose [mg/dL]

700 350 350

0 0 4.4

0, 70, 120, 180, 240 0, 70, 120, 180 0, 70, 120, 180

The time variant effects can be addressed by time-variant models. These models consider the change in time of the system by extended models with time variant parameters. The reader can find information regarding these models in [42–44]. However, a well known practical alternative to deal with the time variations, when possible, is to reduce the acquisition time. This can be done by increasing the minimum excited frequency or by reducing the number of experiments. It is important to note that these considerations are not exclusive of the ORPM excitation, but any other excitation signal in time or frequency domain. High amplitudes could trigger nonlinearities while long acquisition times can result in time variant systems. The proposed model in our methodology is a general rational model since it considers an integer model section and a fractional one, as indicates (5). Observe that this model is non-linear in the parameters and requires initial estimation values for the optimization procedure. The proposed methodology includes a procedure to identify such initial parameters based on the integer model. The general method-procedure is summarized in Table 1 along with the consideration here discussed. 3. Results 3.1. Analysis of the BLA of EIS measurements on glucose solutions In order to shows the ability of the methodology to discriminate glucose concentrations at physiological levels, a group of impedance measurements on saline–glucose and saline–albumin–glucose are elaborated. The measurements are developed in the frequency range from 1 Hz to 30 kHz at a controlled temperature of 37◦ and repeated five times for each glucose concentration. Three base solutions are used. The first one is a NaCl solution at 700 mg/dL. This concentration is near the average physiological value. The second solution has a reduced NaCl concentration (350 mg/dL). This concentration represents the minimum physiological, levels of hyponatremia/hipocloremia. The third solution has 350 mg/dL of NaCl concentration and albumin at 4.4 g/dL. Glucose is added until physiological (normal and pathological) concentrations are reached. Table 2 summarizes the used solutions. The solutions were done using demineralized water from a Millipore MilliQ Element system, sodium-chloride (NaCl) at 99% from VWR

Fig. 2. Magnitude of the Zˆ BLA (jωk ) for the solution 1 in Table 2. The colors represent different glucose concentrations. The level of noise (.) and the level of total standard deviation (+) are approximately to the same level indicating that the data were taken at linear or pseudo-linear state. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

laboratory, albumin from human serum >97% and glucose at 99% from Sigma–Aldrich Laboratories. Fig. 2 presents the magnitude of the BLA of the impedance (Zˆ BLA (jωk )) for the solution group 1 in Table 2 (NaCl solution at 700 mg/dL and different glucose concentrations). This figure shows high impedance values at low frequencies. As the frequency increases, the impedance decreases until around 10 kHz. For frequencies higher than >10 kHz, the impedance becomes nearly horizontal. The figure shows the level of noise (ˆ BLA,n (jωk )) and the level of total standard deviation (ˆ BLA,total (jωk )), defined as the levels of noise plus the level of non-linearities. Both of them are at the same level (around −45 dB below the impedance data). This means that the measurements have been acquired in a nearly linear state. The system is stationary along the acquisition time, so the impedance can be accurately modeled by the BLA. Similar results are present for all the impedance solutions in Table 2. Related results have been reported in [45,46], for a wider frequency band. In [47] the impedance changes have been identified at low frequencies. Along the results there exist some differences. This can be explained by the sensors and/or the electrode configuration. The configuration used in our experiments is a three electrodes configuration. This configuration allows measuring potential changes of the working electrode, independent from changes that may occur at the counter electrode [5,48]. The experiments were done using screen-printed electrodes from DropSens [49]. Counter and working electrode in Platinum and reference electrode in silver.

184

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

followed by a quasi horizontal line at high frequencies (higher than 10 kHz). Based on this, it seems appropriate to choose two different groups of poles and zeros that follow the inequalities (7) and (8). The group that represents the descendant slope will be compressed to a fractional pole, and the group representing the quasi flat section will be compressed to a fractional zero. The poles and zeros fulfilling the inequalities (7) and (8) are highlighted by an ellipse in the figure. Based on the identified chains, the fractional model (5) becomes Zfrac (ω) = K

Fig. 3. Poles and zeros locations of the integer model over the Bode diagram of a group of experimental impedance data.

(jω − zfrac )

ˇ

(jω − pfrac )

˛

(13)

Note that in this model the integer section GR (jω) in (5) is equal to one, since all the poles and zeros have been compressed in the fractional structure (13). The compression equations (9) to (11) generate a first estimation of the parameter values, while the optimization process (12) gives the final values.

3.2. Modeling EIS data by general fractional model In this section the glucose discrimination using fractional models is presented together with the parameters behavior and their relation with the glucose content. 3.2.1. Fractional model identification As has been presented in Section 2.2, the identification based fractional models starts with an integer model estimation in order to obtain a first assessment of the fractional structure and its parameters. This identification procedure generates an integer order model with 7 poles and 7 zeros (7/7). The integer model selection has been chosen based on the Akaike’s information criterion [50], the minimum description length [51], and residues behavior. From the integer model, the poles and zeros present intercalated locations suggesting that the system is prone to be modeled by a fractional model structure, as has been previously illustrated in Section 2.2). Fig. 3 illustrates the impedance response for the solution group 1 in Table 2 and the location of the poles and zeros of the integer model. The behavior of the impedance response can be described as a descendant slope at low frequencies (<10 kHz)

3.2.2. Parameters behavior In order to identify whether the fractional model discriminates the glucose content, the parameters of the model are analyzed regarding the glucose content for each solutions in Table 2. Each parameter is indicated together the standard error bound (±one standard deviation). Figs. 4–6 show the behavior of the pole and zero of the fractional model for the impedance of the solutions in Table 2. The figures represent the root location in rad/s versus the glucose concentration in mg/dL. Fig. 4 shows the fractional pole and the fractional zero behavior of the impedance model for the solution group 1 in Table 2. From the figure, a clear relation between the fractional pole–zero and the glucose concentration is present. The error bound is small enough to identify and discriminate the different glucose concentrations along the solution. Fig. 5 shows the behavior of the fractional pole and zero of the impedance model for the solution group 2 in Table 2. This solution has low concentration of NaCl (lower conductivity that the solution group 1). Nevertheless, even with the change in conductivity the pole and zero location of the fractional model allow discrim-

Fig. 4. Mean parameter value (*) and ±one standard deviation from the fractional model for the solution 1 in Table 2.

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

185

Fig. 5. Mean parameter value (*) and ±one standard deviation from the fractional model for the solution 2 in Table 2.

Fig. 6. Mean parameter value (*) and ±one standard deviation from the fractional model for the solution 3 in Table 2.

inating the different glucose concentrations. Observe that for the pure saline (no-glucose) the pole location do not follow the tendency between the different glucose concentration. However, in the presence of glucose the tendency is clear, with a significant difference as indicates the error bound, specially in the fractional zero. Fig. 6 presents the pole and zero behavior for the fractional model of the impedance for the group of solutions 3 in Table 2. This solution presents low level of NaCl and albumin at average physiological level (note that the albumin concentration is one order of magnitude higher than the NaCl and glucose). From the figure, the pole do not evidence any tendency regarding the glucose content while the error bounds are high enough (overlap each other) that do not allow differentiating the glucose levels. The behavior of the fractional zero, except for the solution without glucose (no-

glucose), presents a correlation with the glucose content with low error bounds that allows discriminating the different glucose levels. Another important observation in Figs. 4–6 is the change in trend of the fractional zero for the solution with albumin (solution 3 – Table 2) regarding to the solutions without albumin (solutions 1 and 2 – Table 2). At this point of our research we do not have a close answer for this behavior however, we are hypothesizing that the change in trend could be given because albumin can bind other molecules [52] and since its concentration is one order of magnitude higher than the other products could affect strongly the physical chemical properties of the analyzed solution. The behavior of the parameters ˛ and ˇ of the fractional model (13) for the impedance data of the solutions in Table 2 are present in Fig. 7. The parameters exhibit small changes at the different glucose concentrations and seem to locate around a specific value for

186

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

Fig. 7. Mean value for the ˇ and ˛ parameters value and ±one standard deviation for the solutions in Table 2. The figure indicates that both parameters present small changes given the glucose content, while the change given the type of solution is relatively high. (For interpretation of the references to color in this figure citation, the reader is referred to the web version of this article.)

the different solutions with glucose. In this sense these parameters could help to identify the type of base solution. For example, for the solution with albumin (red data in Fig. 7) the parameters present the highest values: around 0.82 and 0.95 for ˇ and ˛ respectively. For the solution without albumin and low NaCl content (green data in Fig. 7) the parameters ˇ and ˛ are the lowest ones, around 0.74 and 0.9 respectively. For the solution with the highest concentration of NaCl, the ˇ and ˛ parameters were located at intermediate values around 0.78 and 0.91 respectively. In order to quantify the change in parameters ˇ and ˛ their percentage of variation is calculated (variation of 100% is a change covering the full range 0–1). For the solution 1, the parameter ˇ presents a percentage of variation of 4%, while the parameter ˛ shows a percentage of variation around 2%. For solution 2 the percentage of variation are 2.4% and 1.5% for ˇ and ˛ respectively. For the solution 3 the percentage of variation is 0.77% for ˇ and 0.8% for ˛. 4. Glucose discrimination using circuit and integer models – comparison In this section the impedance modeling is done using circuit and integer models. The results allow comparing the methods showing that even when all the methodologies are sensible to changes in the impedance, given glucose changes, the fractional model allows discriminating the glucose content. A extended explanation of the integer and circuit models can be found in [34,53]. In the current section we limit to a introduction description of the models and the comparison regarding the solutions described in Table 2. 4.1. Circuit models Circuit models is a pragmatic methodology to model and identify changes in the impedance. In order to develop a circuit model it is important to know deeply the processes or phenomena that occur in the analyzed system. When an electrode is in contact with a physiological medium (interstitial fluid, blood, intra/extracellular

Fig. 8. Circuit model – (CPE/Rct ) − Rhf .

fluid, etc.) redistribution of charges is present and an interface is formed [54]. This interface can be modeled by a distributed capacitor known as constant phase element (CPE). The CPE is defined as ZCPE (jω) =

1 Q (jω)



(14)

where 1/Q is a measure of the CPE impedance and has units of F · s−1 , while  is an empirical constant ranging from 0 to 1 [55]. In the interface not only redistribution of charges are present, but a transfer could occurs between the electrode and electrolyte. This process can be modeled by a parallel resistance given the transfer nature of the phenomena Rct [56]. The described interface impedance is located physically in series with the electrolyte, which is represented by a high frequency limiting resistance Rhf [57,58]. Fig. 8 shows the circuit model here used. In the circuit model the behavior of the parameters regarding the glucose concentration is analyzed since they are related with a physical/chemical interpretation. The estimated parameter values of the circuit model along the different group of solutions and glucose concentrations are plotted in Figs. 9 and 10 . Fig. 9 presents the behavior of the Q and  parameters of the CPE defined in Eq. (14). The parameter Q for the saline solutions at 700 mg/dL (solution 1, Table 2) shows increments when the glucose content increases. However, for the other solutions (solution 2 and 3, Table 2), taking into account the error bound, such relation is not well established. Regarding the  parameter, it does not allow dif-

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

187

Fig. 9. Behavior of the parameters Q and  of the CPE element along the different solutions and glucose concentrations.

Fig. 10. Behavior of the parameters Rct and Rhf along the different solutions and glucose concentrations.

ferentiating the glucose concentrations for any of the solutions, but shows a clear difference between the solutions with and without albumin. For the solutions without albumin the  value is located around 0.96, while for the solution without albumin, the parameter value is located around 0.91. Fig. 10 shows the behavior of the parameters Rct and Rhf for the circuit model. The parameter Rct shows tendency for the glucose content in solution 1 (pure saline solution at 700 mg/dL). However, once the level of NaCl is reduced (solution 2) or in the presence of albumin (solution 3) the parameter does not allows discriminating

the glucose content. Regarding the parameter Rhf , the figure shows that the parameter value increases at glucose increments for the solution 1. However, for the other solutions in Table 2 the parameter does not shows a clear tendency and/or the error bound does not allows discriminating between concentrations. From the results it is possible to establish that changes in glucose concentration affect the impedance, but in general, the circuit models do not reveal a clear relations between the parameters behavior and the glucose content in the solutions.

188

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

Fig. 11. Mean zero (top) and pole (bottom) locations of the rational integer model for the group of solutions 1 in Table 2. The colors refer the glucose content. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 12. Mean position plus ±one standard deviation for the zeros nos. 2, 3 and 7 (top) and poles nos. 2, 3, 6 (bottom) from the model of the impedance from the solution 1 (saline solution at 700 mg/dL) in Table 2.

4.2. Integer models In this section a general integer rational model is used. The model identified has been previously introduced in Eq. (6) and repeated here by convenience:

mB (s − zm ) m=1 G(s) = K m A m=1

(s − pm )

(15)

The rational integer model identification leads to a model of order mB = 7 and mA = 7. As example of the behavior of the poles and zeros in function of the glucose content, Fig. 11 shows their mean locations for the solution 1 in Table 2. The colors in the figure

refers the glucose concentration. From the figure, it is possible to identify that the poles and zeros are grouped and their positions are influenced by the glucose. Similar behavior is obtained for the other solutions in Table 2. The poles and zeros of the integer model are analyzed as a function of the glucose in a similar way as was presented for the fractional model. Given the high number of poles and zeros the following convention is used: the couple of pole and zero furthest from the origin is denoted as pole or zero number one (no. 1). The previous one is denoted as pole/zero no. 2, and so on. In this sense the dominant pole and zero, the closest to the origin, are denoted as pole/zero no. 7.

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

189

Fig. 13. Mean position plus ±one standard deviation for the zeros nos. 2, 3 and 7 (top) and poles nos. 2, 3, 6 (bottom) from the model of the impedance from the solution 2 (saline solution at 350 mg/dL) in Table 2.

Fig. 14. Mean position plus ±one standard deviation for the zeros nos. 2, 3 and 7 (top) and poles nos. 2, 3, 6 (bottom) from the model of the impedance from the solution 3 (saline solution at 350 mg/dL and albumin at 4.4 g/dL) in Table 2.

As illustration of the behavior of the poles and zeros of the integer model regarding the glucose content, Fig. 12 presents the mean position together the error bound (±one standard deviation) for the zeros nos. 2, 3 and 7 (top) and poles nos. 2, 3, 6 (bottom) using the impedance data from the solution group 1 in Table 2. From the figure it is possible to note that zeros 3 and 7 (top, Fig. 12) present a slight tendency regarding the glucose. However, once the standard deviation is considered it is not possible to discriminate with confidence the different concentrations. Regarding the poles (bottom, Fig. 12), it seems that they do not show any trend as a function of

the glucose. For the other poles and zeros (not shown in the figure) no clear relation with the glucose content has been identified. If well the zeros 3 and 7 for the solution 1 present a slight correlation with the glucose, such behavior is not present along the other solutions. In general, considering the error bound, the position of the poles and zeros seems unaltered by the glucose content. Figs. 13 and 14 present the behavior of the zeros nos. 2, 3 and 7 (top) and the poles nos. 2, 3, 6 (bottom) of the integer model for the solutions 2 and 3 in Table 2. The figures do not show a readable trend between the pole/zero location and the glucose content. For

190

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

the other poles and zeros, not shown in the figure, similar results were obtained. 5. Conclusion The manuscript addresses the problem of discriminating the glucose at physiological levels using EIS data. The proposed methodology is based in a two steps process: appropriate excitation signal strategy and a fractional model identification. The convenient excitation will allow easily assessing the condition of linearity and time-invariance along the acquisition time. Although the assessment of these conditions seems a basic step, frequently they are taken for granted inducing to errors in the impedance analysis. The fractional model identification is proposed based on the characteristic behavior of the impedance data and the roots location in the integer model in the Laplace domain. The integer model presents high order complexity while its poles and zeros are intercalated making the system prone to be modeled by a fractional model. In this sense, the use of a fractional models is convenient since reduces the total complexity of the model and from here the parameters uncertainty. This allows discriminating the glucose content at physiological levels. The fractional model shows a tight correlation between the glucose content and the fractional zero of the model. The methodology has been analyzed using glucose concentrations at physiological levels in hypo-, normo- and hyper-glycemia. Two different levels of NaCl have been selected in order to obtain different conductivity conditions, and finally, albumin with an order of magnitude higher than the glucose. These solutions if well do not include all blood/plasma properties and characteristics allows presenting a general view the robustness of the methodology here presented. The results obtained with fractional models have been compared with circuit models and rational integer models showing the advantages of the proposed general fractional model identification methodology. From the different modeling methodologies it is clear that glucose, at physiological levels, affects the impedance spectra. However, circuit models and integer rational models were not able to discriminate the glucose changes. The principal difficulty in these models is the high uncertainty introduced for the ordercomplexity of the model, or because the model structure requires accurate knowledge of the processes under analysis. On the other hand, Fractional models present a reduced number of parameters while the initial values of the identification process can be easily identified based on the structure of the integer Laplace model. The result shows the competence of fractional models to identify the glucose changes over impedance measurements at different scenarios: high and low conductivity and in presence of albumin. Based on these results, fractional models present a feasible alternative to interpret the impedance behavior of biological solutions, and a tool for the follow-up of glucose at physiological levels using EIS technology. References [1] K. Tonyushkina, J. Nichols, Glucose meters: a review of technical challenges to obtaining accurate results, J. Diabetes Sci. Technol. 3 (4) (2009) 971–980. [2] T. Koschinsky, L. Heinemann, Sensors for glucose monitoring: technical and clinical aspects, Diabetes/Metab. Res. Rev. 17 (2) (2001) 113–123, http://dx. doi.org/10.1002/dmrr.188. [3] D.C. Klonoff, A review of continuous glucose monitoring technology, Diabetes Technol. Therap. 7 (5) (2005) 770–775, http://dx.doi.org/10.1089/dia.2005.7. 770. [4] O.S. Khalil, Non-invasive glucose measurement technologies an update from 1999 to the dawn of the new millennium, Diabetes Technol. Therap. 6 (5) (2004) 660–697. [5] V. Lvovich, Impedance Spectroscopy: Applications to Electrochemical and Dielectric Phenomena, Wiley, 2012 http://books.google.be/ books?id=CgGqMeQJArkC.

[6] E. Barsoukov, J. Ross Macdonald, Impedance Spectroscopy: Theory, Experiment, and Applications, 2nd ed., Wiley, 2005. [7] S.-M. Park, J.-S. Yoo, Peer reviewed: electrochemical impedance spectroscopy for better electrochemical measurements, Anal. Chem. 75 (21) (2003) 455A–461A, http://dx.doi.org/10.1021/ac0313973. [8] S. Jensen, A. Hauch, P. Hendriksen, M. Mogensen, N. Bonanos, T. Jacobsen, A method to separate process contributions in impedance spectra by variation of test conditions, Electrochem. Soc. J. 154 (12) (2007) B1325–B1330, http:// dx.doi.org/10.1149/1.2790791. [9] F. Groeber, L. Engelhardt, S. Egger, H. Werthmann, M. Monaghan, H. Walles, J. Hansmann, Impedance spectroscopy for the non-destructive evaluation of in vitro epidermal models, Pharmaceut. Res. 32 (5) (2015) 1845–1854, http:// dx.doi.org/10.1007/s11095-014-1580-3. [10] P. Aberg, P. Geladi, I. Nicander, J. Hansson, U. Holmgren, S. Ollmar, Non-invasive and microinvasive electrical impedance spectra of skin cancer, a comparison between two techniques, Skin Res. Technol. 11 (4) (2005) 281–286, http://dx.doi.org/10.1111/j.0909-725X.2005.00125.x. [11] A.R. West, D.C. Sinclair, N. Hirose, Characterization of electrical materials, especially ferroelectrics, by impedance spectroscopy, J. Electroceram. 1 (1) (1997) 65–71, http://dx.doi.org/10.1023/A:1009950415758. [12] O. Olarte Rodriguez, K. Barbé, W. Van Moer, Y. Van Ingelgem, Fractional frequency domain identification of NaCl–glucose solutions at physiological levels, Measurement 50 (2014) 213–221. [13] R. Pintelon, J. Schoukens, System Identification: A Frequency Domain Approach, Wiley and IEEE Press, 2012. [14] A. Elwakil, Fractional-order circuits and systems: an emerging interdisciplinary research area, IEEE Circuits Syst. Mag. 10 (4) (2010) 40–50, http://dx.doi.org/10.1109/MCAS.2010.938637. [15] M. Ortigueira, An introduction to the fractional continuous-time linear systems: the 21st century systems, IEEE Circuits Syst. Mag. 8 (3) (2008) 19–26. [16] M.D. Ortigueira, Fractional Calculus for Scientists and Engineers, Vol. 84 of Lecture Notes in Electrical Engineering, Springer, 2011. [17] J. Sabatier, O.P. Agrawal, J.A.T. Machado, Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering, 1st ed., Springer Publishing Company Incorporated, 2007. [18] K. Barbé, O. Olarte, W.V. Moer, L. Lauwers, Fractional models for modeling complex linear systems under poor frequency resolution measurements, Dig. Signal Process. 23 (4) (2013) 1084–1093, http://dx.doi.org/10.1016/j.dsp. 2013.01.009. [19] K. Barbé, W. Van Moer, G. Nagels, Fractional-order time series models for extracting the haemodynamic response from functional magnetic resonance imaging data, IEEE Trans. Biomed. Eng. 59 (8) (2012) 2264–2272, http://dx. doi.org/10.1109/TBME.2012.2202117. [20] C. Ionescu, J.T. Machado, R.D. Keyser, Fractional-order impulse response of the respiratory system, Comput. Math. Appl. 62 (3) (2011) 845–854 (Special Issue on Advances in Fractional Differential Equations II). [21] J.T. Machado, A.C. Costa, M.D. Quelhas, Fractional dynamics in {DNA}, Commun. Nonlinear Sci. Numer. Simul. 16 (8) (2011) 2963–2969, http://dx. doi.org/10.1016/j.cnsns.2010.11.007. [22] A. Arafa, S. Rida, M. Khalil, Fractional modeling dynamics of HIV and CD4+ T-cells during primary infection, Nonlinear Biomed. Phys. 6 (1) (2012) 1–7, http://dx.doi.org/10.1186/1753-4631-6-1. [23] F.A. Rihan, Numerical modeling of fractional-order biological systems, Abstr. Appl. Anal. 2013 (2013), http://dx.doi.org/10.1155/2013/816803, 11 pp. [24] R.L. Magin, Fractional calculus models of complex dynamics in biological tissues, Comput. Math. Appl. 59 (5) (2010) 1586–1593, http://dx.doi.org/10. 1016/j.camwa.2009.08.039, Fractional Differentiation and Its Applications. [25] M.E. Orazem, T. Bernard, Electrochemical Impedance Spectroscopy, Wiley, ECS The Electrochemical Society Series, 2008. [26] L. Pauwels, W. Simons, A. Hubin, J. Schoukens, R. Pintelon, Contribution to the impedance data analysis of mass-transfer controlled electrochemical systems, in: New Trends in Electrochemical Impedance Spectroscopy (EIS) and Electrochemical Noise Analysis (ENA), 2001, pp. 23–38. [27] T. Freeborn, A survey of fractional-order circuit models for biology and biomedicine, IEEE J. Emerg. Sel. Top. Circuits Syst. 3 (3) (2013) 416–424, http://dx.doi.org/10.1109/JETCAS.2013.2265797. [28] J.T. Machado, V. Kiryakova, F. Mainardi, Recent history of fractional calculus, Commun. Nonlinear Sci. Numer. Simul. 16 (3) (2011) 1140–1153, http://dx. doi.org/10.1016/j.cnsns.2010.05.027. [29] J. Schoukens, R. Pintelon, Y. Rolain, Broadband versus stepped sine FRF measurements, IEEE Trans. Instrum. Meas. 49 (2) (2000) 275–278, http://dx. doi.org/10.1109/19.843063. [30] Y. Van Ingelgem, E. Tourwe, O. Blajiev, R. Pintelon, A. Hubin, Advantages of odd random phase multisine electrochemical impedance measurements, Electroanalysis 21 (6) (2009) 730–739, http://dx.doi.org/10.1002/elan. 200804471. [31] T. D’haene, R. Pintelon, J. Schoukens, E. Van Gheem, Variance analysis of frequency response function measurements using periodic excitations, IEEE Trans. Instrum. Meas. 54 (4) (2005) 1452–1456, http://dx.doi.org/10.1109/ TIM.2005.851075. [32] E. Van Gheem, A New Methodology For Electrochemical Impedance Spectroscopy In The Presence Of Non-Linear Distortions And Non-Stationary Behaviour: Application to Pitting Corrosion of Aluminium (Ph.D. thesis), Vrije Universiteit Brussel VUB, 2005. [33] E. Van Gheem, R. Pintelon, J. Vereecken, J. Schoukens, A. Hubin, P. Verboven, O. Blajiev, Electrochemical impedance spectroscopy in the presence of

O. Olarte, K. Barbé / Biomedical Signal Processing and Control 40 (2018) 180–191

[34]

[35]

[36]

[37]

[38]

[39] [40]

[41]

[42]

[43]

[44]

non-linear distortions and non-stationary behaviour. Part I. Theory and validation, Electrochim. Acta 49 (26) (2004) 4753–4762, http://dx.doi.org/10. 1016/j.electacta.2004.05.039. O. Olarte Rodriguez, K. Barbé, W. Van Moer, Y. Van Ingelgem, A. Hubin, Measurement and characterization of glucose in NaCl aqueous solutions by electrochemical impedance spectroscopy, Biomed. Signal Process. Control 14 (2014) 9–18. J. Schoukens, R. Pintelon, T. Dobrowiecki, Y. Rolain, Identification of linear systems with nonlinear distortions, Automatica 41 (3) (2005) 491–504, http://dx.doi.org/10.1016/j.automatica.2004.10.004 (Data-Based Modelling and System Identification). J. Schoukens, R. Pintelon, Y. Rolain, T. Dobrowiecki, Frequency response function measurements in the presence of nonlinear distortions, Automatica 37 (6) (2001) 939–946, http://dx.doi.org/10.1016/S0005-1098(01)00037-1. R. Mansouri, M. Bettayeb, S. Djennoune, Approximation of high order integer systems by fractional order reduced-parameters models, Math. Comput. Modell. 51 (12) (2010) 53–62, http://dx.doi.org/10.1016/j.mcm.2009.07.018. A. Charef, H.H. Sun, Y.-Y. Tsao, B. Onaral, Fractal system as represented by singularity function, IEEE Trans. Autom. Control 37 (9) (1992) 1465–1470, http://dx.doi.org/10.1109/9.159595. K. Worden, R. Tomlinson, Nonlinearity in Structural Dynamics. Detection, Identification and Modelling, IOP Publishing Ltd., London, 2001. G. Kerschen, K. Worden, A. Vakakis, J. Golinval, Past, present and future of nonlinear system identification in structural dynamics, Mech. Syst. Signal Process. 20 (3) (2006) 505–592, http://dx.doi.org/10.1016/j.ymssp.2005.04. 008. J.P. Noël, G. Kerschen, Nonlinear system identification in structural dynamics: 10 more years of progress, Mech. Syst. Signal Process. 83 (2017) 2–35, http:// dx.doi.org/10.1016/j.ymssp.2016.07.020. J. Lataire, E. Louarroudi, R. Pintelon, Detecting a time-varying behavior in frequency response function measurements, IEEE Trans. Instrum. Meas. 61 (8) (2012) 2132–2143, http://dx.doi.org/10.1109/TIM.2012.2197072. R. Pintelon, E. Louarroudi, J. Lataire, Detection and quantification of the influence of time variation in frequency response function measurements using arbitrary excitations, IEEE Trans. Instrum. Meas. 61 (12) (2012) 3387–3395, http://dx.doi.org/10.1109/TIM.2012.2210327. K. Liu, L. Deng, Experimental verification of an algorithm for identification of linear time-varying systems, J. Sound Vibr. 279 (3) (2005) 1170–1180, http:// dx.doi.org/10.1016/j.jsv.2004.01.040.

191

[45] J.H. Park, C.S. Kim, B.C. Choi, K.Y. Ham, The correlation of the complex dielectric constant and blood glucose at low frequency, Biosens. Bioelectron. 19 (4) (2003) 321–324, http://dx.doi.org/10.1016/S0956-5663(03)00188-X. [46] A. Tura, S. Sbrignadello, S. Barison, S. Conti, G. Pacini, Impedance spectroscopy of solutions at physiological glucose concentrations, Biophys. Chem. 129 (2–3) (2007) 235–241, http://dx.doi.org/10.1016/j.bpc.2007.06.001. [47] R.K. Shervedani, A.H. Mehrjardi, N. Zamiri, A novel method for glucose determination based on electrochemical impedance spectroscopy using glucose oxidase self-assembled biosensor, Bioelectrochemistry 69 (2) (2006) 201–208, http://dx.doi.org/10.1016/j.bioelechem.2006.01.003. [48] T.N. Gamry Instruments, Two-, Three-, and Four-Electrode Experiments, Tech. Rep., Gamry Instruments. http://www.gamry.com/application-notes/allnotes/. [49] DropSens, Screen-Printed Electrodes, Tech. Rep. (accessed 28.01.14). [50] H. Akaike, New look at statistical-model identification, IEEE Trans. Autom. Control AC19 (6) (1974) 716–723. [51] J. Rissanen, Modeling by shortest data description, Automatica 14 (5) (1978) 465–471. [52] Sigma-Aldrich, Product Information, Albumin, Human, Tech. Rep., Sigma Aldrich. [53] O. Olarte Rodriguez, K. Barbé, W. Van Moer, Y. Van Ingelgem, Glucose Characterization Based on Electrochemical Impedance Spectroscopy, IEEE Instrumentation and Measurement Society, 2014, pp. 833–837. [54] E. McAdams, Encyclopedia of Medical Devices and Instrumentation, 2nd ed., John Wiley and Sons Inc., 2006. [55] R.L. Hurt, J. Macdonald, Distributed circuit elements in impedance spectroscopy: a unified treatment of conductive and dielectric systems, Solid State Ionics 20 (2) (1986) 111–124, http://dx.doi.org/10.1016/01672738(86)90018-4. [56] D.R. Merrill, M. Bikson, J.G. Jefferys, Electrical stimulation of excitable tissue: design of efficacious and safe protocols, J. Neurosci. Methods 141 (2) (2005) 171–198, http://dx.doi.org/10.1016/j.jneumeth.2004.10.020. [57] V. Lvovich, S. Srikanthan, R.L. Silverstein, A novel broadband impedance method for detection of cell-derived microparticles, Biosens. Bioelectron. 26 (2) (2010) 444–451, http://dx.doi.org/10.1016/j.bios.2010.07.094. [58] J. Macdonald, Impedance spectroscopy, Ann. Biomed. Eng. 20 (3) (1992) 289–305, http://dx.doi.org/10.1007/BF02368532.