Investigating the electrode-electrolyte interface modelling in cochlear implants

Objective. Proposing a good electrode-electrolyte interface (EEI) model and properly identifying relevant parameters may help designing safer and more optimized auditory nerve fiber stimulation and recording in cochlear implants (CI). However, in literature, EEI model parameter values exhibit large variability. We aim to explain some root causes of this variability using the Cole model and its simpler form, the Basic RC model. Approach. We use temporal and spectral methods and fit the models to stimulation pulse voltage response (SPVR) and electrochemical impedance spectroscopy (EIS) data. Main Results. Temporal fittings show that there are multiple sets of model parameters that provide a good fit to the SPVR data. Therefore, small methodological differences in literature may result in different model fits. While these models share similar characteristics at high frequencies >500 Hz, the SPVR fitting is blind to low frequencies, thus it cannot correctly estimate the Faradaic resistor. Similarly, the polarization capacitor and its fractional order are not estimated robustly (capacitor variations in the nano- to micro-farad range) due to limited observation of mid-range frequencies. EIS provides a good model fit down to ∼3Hz, and thus robust estimation for the polarization capacitor. At lower frequencies charge mechanisms may modify the EEI, requiring multi-compartment Cole model fitting to EIS to improve the estimation of Faradaic characteristics. Our EIS data measurements down to 0.05Hz show that a two-compartment Cole model is sufficient to explain the data. Significance. Our study describes the scope and limitation of SPVR and EIS fitting methods, by which literature variability is explained among CI EEI models. The estimation of mid-to-low-frequency characteristics of the CI EEI is not in the scope of the SPVR method. EIS provides a better fit; however, its results should not be extrapolated to unobserved frequencies where new charge transfer mechanisms may emerge at the EEI.


Background and the study purpose
Studying the characteristics of the electrode-tissue and the electrode-electrolyte interface (EEI) by means of mathematical models helps to improve our understanding in several research areas, including neuromodulation-related technologies such as cochlear implants (CI). Specific to CIs, these models aim to enhance our knowledge on stimulation current spread, neuronal excitability along the cochlea, and auditory nerve response [1,2]. A reliable interface model also supports designing efficacious and safe stimulation and acquisition protocols [3]. Widely accepted models for the electrode-tissue interface or EEI include an ohmic resistor R S , in series with a complex-value polarization impedance Z P . The impedance Z P is usually composed of a Faradaic process polarization resistor, R P , in parallel with a double-layer perfect capacitor C P or with a constant phase element (CPE) characterized by an imperfect (pseudo) capacitor with fractional order α [4]. These models are graphically depicted in figure 1, and are called hereafter the Basic and the Cole model, respectively. The electrode roughness and its fractal dimension, inhomogeneous reaction rates on the electrode surface, crystal orientation, and resistance distribution in an oxide layer at the surface are some of the physical explanations put forth for a CPE [5,6]. It has been shown that a capacitive CPE is modeled with a timevarying inductance in series with a resistor [7]. This resistor and its corresponding energy dissipation explains imperfectness of CPE as a capacitor.
A great deal of cochlear implant-related literature has examined and proposed electrode-tissue interface and EEI models based on platinum electrodes tested in vivo and in vitro using an electrolyte mimicking human perilymph [8][9][10][11][12][13][14][15][16]. These studies characterize and model the EEI using both electrochemical impedance spectroscopy (EIS) and stimulation pulse voltage response (SPVR) data. A wide range of parameter values are reported, with variability over three orders of magnitude for C P and R P (e.g., R P is equal to a few kΩ in [10], to a few hundreds of kΩ in [16], and to >10 12 kΩ in [15]).
In this work, we examined the state-of-the-art to understand literature discrepancies in terms of models and parameter values. We investigated the large divergences and their root causes. To achieve this goal, we carried out a well-controlled in vitro data recording using CI electrodes and investigated the Basic and Cole models fittings (i.e., the two widely accepted models for CI EEI modeling) in both the time and frequency domains.
1.2. State-of-the-art: inconsistencies EEI models have been developed and characterized in literature using both time and frequency domains with SPVR and EIS techniques. We summarize the details and differences found in literature and provide an overview of the varying reported values.

Temporal EEI modelling
Lim et al [8] modelled the cochlear implant EEI using the Basic model by means of in vivo measurements using biphasic current pulse of 172.8μs/phase. They showed that R S , R P , and C P are approximately equal to 3kΩ, 5kΩ and 10nF when the applied current to the electrode (width: 0.3mm, ⌀0.4 or ⌀0.6 mm) was 0.47 mA (≈1 mA mm −2 ). Vanpoucke et al [9] simulated a CI electrode array (HiFocus electrode array with rectangular platinum contacts measuring 0.4×0.5 mm 2 spaced on a 1.1 mm pitch) by placing 16 Basic models in a resistive ladder network that mimics the intracochlear structure. The model parameters were identified in vivo using 66 μs/phase pulse stimulation and changing the current levels from 40 to 400μA. The parameter values slightly varied with patient, stimulating current value, and electrode number. Generally, the bulk resistor was in the range of 2-4kΩ, C P was approximately 5nF, and R P was on the order of 5-7kΩ.
Tykocinski et al [10] performed a clinical study and modeled the EEI for both the conventional fullbanded straight and Contour versions of the Nucleus-24 electrode array (Cochlear Ltd.), with electrodes sized 0.48 mm 2 and 0.3 mm 2 , respectively. They used 76μA, 25μs/phase monopolar biphasic stimulation. The mean values of R P ,C P were about 2kΩ,9nF and 4kΩ,5nF for the conventional and the contour electrodes, respectively. The R S was ∼3.5 kΩ. is obtained from the Cole model when the CPE component is substituted by a perfect capacitor C P having a non-fractional order (i.e., α = 1). A two compartment Cole model (right) is obtained by putting two Cole models in series. The pseudo-capacitor coefficient C P in the Cole model has the units of Ω −1 (jω) −α and its value is expressed with F·sec α−1 (where F stands for Farads). For simplicity we use the abbreviation F * to refer to this unit. Di Lella et al [13] identified the parameters of the Basic model using clinical CI telemetry voltage acquired in different grounding configuration modes (Cochlear Nucleus Contour Advance electrode). The mean values of R S , R P , and C P were approximately equal to 3.3kΩ, 7kΩ, and 6nF, respectively. This method assumed that no current flows through the Faradaic resistor shortly after the current step is applied. Although this assumption facilitates approximate calculation of C P , the results may be sensitive to the sampling rate and noise level of the recorded signals. One year after this study, Di Lella et al [17] aimed to validate their proposed method using emulated implant loads. The simulated R S and C P values could be estimated well, however the value of R P could not due to negligible variation of the Z P during the short application time of the current source (<60 μs). Therefore, although R P could be potentially used as an assessment parameter for CI health, the authors concluded that the most important subcomponents to be considered for future analysis are R S and C P .

Spectral EEI modelling
To our knowledge, the research conducted by Mesnildrey et al [15] was the first effort to fit a cochlear implant EEI (Advance Bionics ® HiFocus 1J electrode array, consists of 16 planar contacts of 0.4×0.5 mm 2 ) with a Cole model using temporal and spectral methods. Both in vitro and clinical experiments were performed, and it was found that cochlear tissue can be well represented by a pure resistive impedance. They showed that onset, offset, and phase reversals of temporal responses are smoothed exponentially due to a parasitic capacitor (∼0.34 nF) that is formed between electrode wires and the power line. This finding suggests that environmental parasitic factors may potentially affect the results of fitting methods that are based on morphological analysis of the responses at onset or offset transient (as used in [13]).
Mesnildrey et al removed the Faradaic resistor R P from the model because their in vitro data suggested a very big R P value (>10 15 Ω). This estimation is not comparable with any preceding EEI modeling study conducted in CI, although it somehow is in agreement with the R P = 38.5 MΩ reported by Zhou et al [14] in deep brain stimulation using ⌀1mm Pt electrodes. Concerning the CPE parameters, the mean values of α and C P were estimated to be 0.64 and 481·10 −9 Ω -1 (jω) -α in the time domain, and ∼0.75 and ∼200·10 −9 Ω -1 (jω) -α (deduced from Figure 10) in the frequency domain, respectively. The pseudo-capacitor C P has the units of Ω -1 (jω) -α and its value is expressed with F·sec α-1 . Hereafter, for simplicity we use the abbreviation F * instead of F·sec α-1 .
Jiang et al [16] aimed to characterize the human cadaveric cochlea complex transimpedance using EIS data from 10 Hz to 100 kHz (Advance Bionics ® HiFocus 1J electrode array). They fit the data by means of a two-compartment Cole model to integrate spread-independent (i.e., EEI) and spread-induced (cochlea tissue) impedances. Regarding their EEI model, the values of R P , C P and α were estimated to be about 120kΩ, 2μF * , and 0.8, respectively. The spreadinduced compartment (i.e., cochlea) showed a resistive characteristic below ∼10 kHz (a pure resistive R S as the bulk) and turned into a complex impedance with increasing frequency due to a CPE component (C P = 1-2 nF * , α = 0.8) used to express the cochlea tissue and/or parasitic effects at high frequencies.
The measurements were collected in an artificial perilymph (without human serum albumin) to simulate the human perilymph. The solution was composed of (in mM) NaCl (125), KCl (3.5), NaHCO 3 (25), CaCl 2 (1.3), MgCl 2 (1.2), and NaH 2 PO 4 (0.75) [18]. The solution was deoxygenated before and during the experiments using N 2 bubbling. The reference and the working electrodes were put as close as possible in the electrolyte.

Chronopotentiometry
Chronopotentiometry was performed to record the SPVR data following a 20sec zero-current injection for the stabilization of the recording. Chronopotentiometry stimulation pulses were symmetric, biphasic, current-controlled, anodic-first pulses with no interphase interval and with a phase duration of 180 μs. For each electrode, the recording was repeated for a set of current amplitude values equal to {100, 250, 500, 750, 1000, 1250, 1500, 1750, 2000}μA. Voltage responses were acquired at a sampling frequency of 500 kHz.

Electrochemical Impedance Spectroscopy
The EIS data (magnitude and phase) was measured at logarithmic-spanned frequencies ranging from 0.05Hz to 1.0MHz (9 samples per decade, giving 67 data points in total). Each measurement was repeated three times and the averaged data was used for an improved signal to noise ratio. This improvement is mainly important for low frequencies below ∼1Hz, where oxidation/reduction reactions may have a larger effect.

Electrode-electrolyte interface model fitting
In this study we model the CI EEI using the Cole model and its simpler form the Basic RC model where the fractional order is equal to 1.0 (see figure 1). The parameters of the EEI model were fit (i.e., estimated) in time and frequency domains using chronopotentiometry (i.e., SPVR) and EIS data, respectively.

SPVR Data Fitting
We used the SPVR data for temporal fitting of the Cole (or Basic) model. We first calculated the theoretical frequency response of the model, then used a numerical method to calculate the voltage response of the model to stimulating current pulses. This process was integrated as a loop in a curve fitting method to find the model parameters that provide the best fit.
The impedance of a generic multi-compartment Cole model (without the effect of bulk resistor) is expressed as follows in the frequency domain: where N is the number of compartments used in the Cole model (e.g., N=2 in the double Cole model in figure 1), and =j 1. We employed the method of linear simulation of fractional-order dynamic systems [19] to approximate the response of the fractional-order Fourier transform ofˆ( ) Z f to input pulse currents. This calculation was conducted by means of the open source Matlab toolbox FOMCON [20]. The model parameters were then estimated by employing a recursive fitting method that minimizes the distance between the model output and the measured SPVR data using the cost function E t expressed as follows: where v model and v meas are the modeled and measured SPVR data, and T ∊{45,90,180}ms is the pulse duration that contributes to the fitting method. In the temporal fitting, the effect of bulk resistor is eliminated (i.e., R S = 0) by removing the abrupt jumps in the measured SPVR data at stimulation pulse onset times; therefore the fit is limited to double layer parameters C P , R P , and α. The double layer or diffuse layer refers to two parallel layers of charge surrounding the electrode, one at the surface of the electrode and one inside the electrolyte with free ions.
In the fitting procedures used in this study, one (or more) double layer parameters may be set to constant value(s). In summary, the following situations may occur: 1. Calculation of E t when all double layer parameters are fixed: This approach is used to obtain 2D (or 3D) topography maps of E t as a function of C P , R P , and α (see figure 2(a)). The maps show which combination of double layer parameters can potentially minimize the cost function and how sensitive they are in generating global and local minima.
2. Fitting C P and R P at fixed α values: This fitting method is used to study the relation between α, R P , and C P and the resulting simulated SPVR data. α is fixed in a range (e.g., 0.5 to 1 with the step size of 0.01) and R P and C P are estimated at different α values (see figure 2(d) and figure 3). Each fitting is executed with 12 different initial conditions obtained from intersecting C P,init = {40,400,4000}nF * and R P,init = {1,10,100,1000}kΩ to reduce the risk of missing the global minimum.
3. All double layer parameters may vary during the fitting process: this approach is used as a generic fitting method. Analogous to item b, the fitting is initiated at 12 different initial conditions. These conditions are intersected with α init = {0.6, 1} to generate 24 total initial conditions.

EIS data fitting
, min max to minimize the cost function E , f expressed as follows: ? ? 3 A normalized form of magnitude difference, as well as phase difference, between the measured and modeled EIS contributes to the cost function E .
f The impedance magnitude is expressed logarithmically so that the cost is not biased with high-impedance lowfrequency terms. The log function acts as a dynamic range compressor function. The final impacts of magnitude and phase on the overall cost function are controlled by coefficients r 1 and r . 2 In this study, r 1 and r 2 were set empirically to 0.75 and 0.25, respectively. Putting more weight on the compressed term facilitates the fit and provides a higher guarantee that | ( )| Z f and |ˆ( )| Z f do not greatly deviate due to a plausible extra compressing effect of the log function.

Results
We firstly report the results of temporal fitting of the Basic and Cole models using SPVR data (section 3.1). A single electrode is first used to explain the details of the concepts and outcomes (as an example). The main results are then given using 10 electrodes. We follow a similar approach for spectral fitting of the models using EIS data (section 3.3). However, before presenting this data, we make a link between these two domains by comparing a measured EIS data to some temporally fitted models (section 3.2). This analysis enabled our determination of the scope and power of using temporally fit models to correctly express characteristics of the EEI.

Temporal fitting
In the temporal fitting we limit the parameters to those of the double layer (i.e., C P , R P , and α). The aim is to investigate how sensitive/robust/reproducible the model output is in response to changes in double layer parameter values and to SPVR data. To do so, we analyze the distribution of low value E t 's when expressed by different sets of double layer parameters. A focal (or spread) distribution indicates that the Cole model is (or is not) sensitive to its parameter variations. Having this information leads us to some of the possible root causes of literature discrepancies.

Single electrode: a detailed case study
In this section the results of the temporal fit are reported for a single electrode.

Topology of the fitting solutions
In this sub-section, we follow the fitting procedure described in section 2.2.1. item a. Figure 2(a) illustrates several 2D color-coded slices of a 3D E t map calculated at different α values for a single electrode. Each map includes a wide range of C P and R P values from 10nF * to 40μF * and from 0.5kΩ to 1TΩ, respectively. For each map, and thus for each given α, there are a horizontal and a vertical area where the cost drops drastically. The cost especially drops along the vertical line, where C P is quasi constant, indicating that the cost function is not sensitive to R P . Increasing the value of α shifts the vertical line to the left, indicating that C P and α are inversely proportional.
When the value of α is big (e.g., α0.9) there is a distinct small area on the map at the intersection of the two lines where the cost is minimal. In this area, the values of C P and R P are relatively low (∼70nF * , ∼2kΩ). When α slightly decreases and takes a value around 0.6-0.7, the spot area moves to the right and changes its shape to an extended area along the vertical line (C P @1000-3000nF * ). This indicates that a Cole model with α@0.65 is almost insensitive to R P . In other words, the cost function E t does not receive enough information from R P to generate a localized minimum area.
Further reduction of the fractional order (e.g., α0.5) increases the cost along the vertical area, and there is a very small gradient toward lowering the cost when R P ¥.
In order to have a better insight on the impact of model parameters on the fitting cost, we reduce the 3D map (along the C P axis) to obtain a 2D map shown in figure 2(b). To do so, the vertical areas from each 2D slice are extracted (i.e., changing 2D to 1D by taking the minimum value of the cost across the C P axis). It is worth mentioning that this dimension reduction does not remove valuable information from the map because C P is quasi constant along the vertical minimum areas, meaning that α and C P are two highly correlated parameters.
The dark blue area in figure 2(b) corresponds to the α and R P values for which the model precisely simulates the SPVR data, and the cost is less than ∼1%. This area is well spread across α (0.6 α 0.9) and R P (3kΩ) axes, meaning that the model is not very sensitive to its parameter values. The sensitivity to R P increases with α. When α = 1.0, a slight deviation of R P from its optimal value around 2kΩ is enough that the cost greatly increases. The least sensitivity is obtained for α@0.63, where the value of R P may change easily from ∼30kΩ to >10 5 kΩ without a significant impact on the cost.
As the model sensitivity to R P is decreased, we expect some shallow local minima to exist in the reduced 2D map of figure 2(b). This is better shown in figure 2(c) where the cost is plotted along the R P axis for various α values. The shallow valleys on the graph corresponding to α = 0.60 and α = 0.65 (see arrows) show a few of these shallow local minima. Figure 2(c) also shows that the global minimum, corresponding to α=0.75, is just 0.6% less than the minimum costs corresponding to α=0.65 and α=0.85. In other words, changing α by 0.1 does not increase the cost E t drastically; therefore, the cost function is not sensitive to the fractional order around its global minimum value. Consequently, slight modifications on the measured SPVR may change the location of this global minimum along the α axis from one fit to another.

The sensitivity and impulse responses
In this sub-section, we follow the fitting procedure described in section 2.2.1 item b. Results indicate that there is a negative quasi-linear relation between α and log(C P ) (figure 2(d)). C P increases from 55nF to 3500nF * when α decreases from 1.0 to 0.6. R P increases at a faster rate with its log value increasing exponentially as α decreases. It is about 2kΩ at α=1, and ∼1TΩ at α=0.6. In figure 2(d) we have highlighted 7 different sets of model parameters for which we compare their simulated SPVRs to the measured SPVR in figure 2(e).
Although these models are quite different in terms of parameter values, their costs are below 2%. Consequently, they are not easily distinguishable from each other and from the measured SPVR, unless the differences to the measured SPVR data is depicted (figure 2(f)). The differences are maximal for the Basic model (#7) and the Cole model corresponding to α = 0.57 (#1). These differences have opposite polarities indicating that the fractional order is over-and under-estimated respectively in these two models. As a result, the curvature of the simulated SPVR by model #7 (#1) is more (less) than the curvature of the measured SPVR.
In models #2, #3, #4, changing the value of R P from a few kΩ to several MΩ at α@0.63 does not significantly change the differences between SPVRs. This is a critical point at which the Cole model has the least sensitivity to R P . Lastly, increasing the α to 0.73 (model #5 where the global minimum is located for this electrode) and to 0.8 (model #6) does not crucially change the differences between SPVRs.

Multi electrode study
In this section the main results of the temporal fit are reported for 10 electrodes.

Sensitivity of the Cole model (cost function) to double layer parameters
We followed the fitting procedure described in section 2.2.1. item b and obtained the relations between α, R P , and C P on 10 electrodes ( figure 3(a)). The results are consistent and similar to the data presented in figure 2(d) (i.e., single electrode data). Log(C P ) and α are related quasi-linearly. C P increases from 66±10nF to 900±140nF * and then to 6700±1000nF * when α decreases from 1.0 to 0.75 and 0.55, respectively. R P increases at a faster rate with its log value increasing exponentially as α decreases. It is equal to 1.7±0.3kΩ, 4.1±0.8kΩ and >1TΩ at the abovementioned αvalues, respectively.
The critical points at which R P increases are around α = 0.63. The cost remains below ∼2% when α varies between 0.6 and 1.0, and minimum values are mostly located between α = 0.7 and α = 0.8 (centered around α=0.75). Modification of α in this range may significantly change the value of C P from ∼500nF * to ∼1500nF * . Figure 3(b) illustrates the boxplots of the parameter values where the cost is minimized for each electrode.
3.1.2.2. The effect of stimulation amplitude and duration on the fit Applying a stimulation current pulse to an electrode at the EEI results in a transient potential at the electrode, required to achieve the defined current level. The transient voltage response is known to be non-linear at the abrupt application of the current pulse [21]. The degree of nonlinearity during the onset of the applied pulse determines if different Cole models are needed to express the system properly. In this section, we study to what extent the parameters of the Cole model are modified by changes in the current amplitude and duration. The answer determines if a generic Cole model can be used in cochlear implant studies. We used current amplitudes ranging from 100 to 2000μA (i.e., ∼0.22 to ∼4.4 mA mm −2 ) and fit Cole models for each amplitude while setting T (i.e., the pulse duration that contributes to the fitting) to 45, 90 and 180μs. Figure 4(a) depicts the modification of the Cole and the Basic models parameters with amplitude and duration. We use the fitting method described in section 2.2.1. item c. The cost for the Cole fit is 0.4%, indicating good fits. The value of α increases with current amplitude for all durations T. It is about 0.78 at 100μA and increases to about 0.9 for big currents 1mA. The value of R P decreases exponentially with current amplitude and duration. When increasing the current from 100μA to 2000μA, R P decreases from 3.1, 5.2 and 7.7kΩ to 0.31, 0.32, and 0.32kΩ, for T =45, 90 and 180μs, respectively. C P does not show a clear increase or decrease tendency. It takes values ranging from ∼0.3μF * to ∼1.5μF * depending on the best determined value of α by the fitting process (for each electrode at a given stimulation amplitude and duration). The standard deviation of α is around 0.1. The model is not so sensitive to the value of α, and this standard deviation is large enough for the estimation of C P with a big standard deviation of ∼0.5μF * (Note: there is a quasi-linear relation between log(C P ) and α for a given current amplitude and duration).
It is noteworthy to mention that the fact there is not a clear decrease in C P despite an increase α with the current amplitude suggests that C P should increase if the value of α is kept constant (like in the Basic model). We verified this by fitting the Basic model to the same data (see figure 4(a) light color dotted lines). The fits are not as good as those of the Cole model. However, the costs remain 2%. As α is always equal to 1, the standard deviations of C P and R P are lower in the Basic model. C P increases linearly with current amplitude and duration, whereas R P decreases exponentially.
The Basic model reveals that when increasing the current amplitude and shortening the time duration, the time constant of the double layer decreases. At 250 μA, when T is set to 45, 90 and 180μs, the time constant (τ = C P ×R P ) is equal to 49±4, 75±5, and 118 ±7μs, respectively. At 1500μA these time constants decrease to 39±4, 47±6, 54±9 μs, respectively. Figure 4(b) illustrates the mean (solid line) and standard deviation (patch area) of the measured and fitted SPVRs using Cole and Basic models fits at two different current amplitudes (250 and 1500μA) and for three different pulse durations T (45, 90, and 180μs). One cannot easily distinguish between the original and fit data (even in the Basic model fit) unless the difference between the two is illustrated. Compared to the Cole model, the Basic model deviates more from the original data in low amplitude (where the best α is around 0.78) than in high amplitude (where the best α is around 0.9). Increasing T increases the deviations of all fit data from the original SPVR data: an indication that not a fully consistent set of model parameters and time constants (especially in the Basic model) can be assigned to different segments (i.e., initial, mid, and end segments) of SPVR data segments.

From time to frequency domain
In this section, we investigate the reasons behind small changes in the simulated SPVR, despite large modifications in model parameter values. To achieve this, impedance spectra (i.e., frequency responses) of the 7 highlighted models in section 3.1.1.2. are calculated and superimposed on the EIS data measured from the same electrode (figure 5). A common characteristic in these models is their impedance spectra similarities at frequencies 500Hz. In this frequency range, EIS data matches the models well. At lower frequencies, the impedance spectra diverge and, in some cases, deviate greatly from the EIS. This indicates that the low-frequency characteristics of the EEI is not considered in the temporal model fitting. In other words, the temporal fitting is blind (semi-blind) to low (medium) frequency characteristics of the EEI. This observation reveals that the estimation of the R P parameter, which is the most important parameter for determining the low frequency characteristics of the Cole/Basic model, is erroneous in the temporal fitting method.
The α parameter, on the other hand, mainly affects the gradient (i.e., tangent of the slope) of the impedance spectrum magnitude at mid-range frequencies.
The spectrum of the Basic model (model #7,α = 1) has the steepest gradient among the Cole models and deviates faster from the EIS data. This has an impact on the estimation of R P and yields in an underestimate of the value (e.g.,2kΩ) compared to what is observed in the EIS data (>1MΩ). Underestimating R P compensates deviation from the EIS at mid-range frequencies around 500 Hz but at the expense of larger deviations for lower frequencies. However, the temporal fitting method is blind to low frequencies and does not increase the cost function.
If the slope of the impedance spectrum of a model and the EIS data are almost in parallel (like in models #2, #3 and #4 at α ∼= 0.63), the upper knee-point on the spectrum can be formed at any frequency range blind to the cost function. Therefore, the value of R P that mostly controls the location of the upper kneepoint can change from 0.8GΩ to35kΩ in models #2 to #4 without any significant impact on the simulated SPVR data (see figures 2(b)-(f)).
If the gradient of the impedance spectrum of a model is smaller than the EIS (like in model #1, α = 0.6), the effect of underestimation of α (i.e., lower gradient) cannot be corrected by R P . This is the α threshold below which R P ¥.
Finally, models #5 and #6 have slightly steeper gradients than the EIS. As a result, their R P values are underestimated (5.9kΩ and 3.8kΩ, respectively) compared to what the EIS data suggests.

Spectral fitting
In the previous section, we showed that the temporal fitting method is blind to low-frequency characteristic of the EEI and mainly explains its high-frequency behavior. In this section, we aim to investigate models that better fit the EIS data and investigate their advantages and limitations. Firstly, using a single electrode, we study in detail the effect of EIS bandwidth on the model fitting outcomes, followed by the investigation on the performance of a two-compartment Cole model in explaining very low frequency characteristics of the EEI (3Hz). The results are then provided for all electrodes.

Single electrode: a detailed case study
In this section the results of the spectral fit are reported for a single electrode.

EIS and the single-compartment Cole model
Before any model fitting is carried out using EIS, important features existing in this data should be well investigated. A good model fitting is the one that can reproduce these features properly. Figure 6 illustrates the single-electrode EIS data that was used previously in section 3.2. This data reveals a resistive (ohmic) characteristic for the EEI at high frequencies (>100 kHz) where the impedance spectrum magnitude starts to saturate at about 500Ω and the phase approaches 0 Rad. EIS mid-range frequencies (∼100Hz) are characterized by a linear slope (indicated as 'Linslope1') of the impedance magnitude and a large negative impedance phase (i.e., capacitive). Mid-range frequencies are basically characterized by the CPE element, while R P and R S have little influence. The slope  Table 1 shows the parameter values of the models. and the horizontal position of the linear region are mostly determined by α and C P , respectively. Below this frequency range, the Cole model expects the EIS magnitude to saturate (i.e., at the R P value) and the phase to increase toward 0 Rad. However, this is not the characteristic of the EIS data presented in figure 6. Below ∼3Hz, the phase returns toward more negative values (creating a bump), and the magnitude deviates to another linear slope (indicated as 'Lin slope 2'). The two slopes are separated by the frequency interval of the quasi-saturation process. These EIS features indicate that a new mechanism is emerged at low frequencies. If fitting a single-compartment Cole model to the EIS is desired, this low-frequency mechanism may be an issue. When very low frequencies (below ∼3Hz) are included in the fitting, they are not perfectly fit and may bias the results. On the other hand, if this mechanism is excluded, the fit model cannot extrapolate its estimations to low frequency terms that have not been observed. This predicament is illustrated in figure 6. The first 5 models share the same single-compartment Cole model architecture. However, they are different in the EIS lower frequency band (f min ) that each model observes in the fitting process: f min ∊{0.05, 0.1, 1, 10, 100}Hz for model #i to #v, respectively. The upper bound is kept constant at 1MHz for all models.
When a model includes low frequencies such that it fully observes the bump and the second linear slope of the EIS (e.g., model #i, #ii), its prediction matches the EIS in a wide frequency range, however at the expense of increased deviations below ∼10kHz. The α parameter is underestimated (see table 1) because a single impedance magnitude slope must be aligned to two linear slopes (i.e., 'Lin slope 1' and 'Lin slope 2') that are separated by a quasi-saturated area. Under these conditions, reduction of the slope helps for a better fit.
When a model does not observe the low-frequency mechanism (e.g., model #iv), the fit is good in the mid-and high-frequency ranges. 'Lin slope 1' and the negative phase 1 are correctly reproduced, however erroneous estimations are obtained at low frequencies because the model wrongly extrapolates (dashed lines in the figure) from the partial observation of the EIS data. Model #iii partially observes the low frequency mechanism (i.e., the bump) and provides a fit that is a tradeoff between those of models #iv and #ii.
A further increase of f min to about 100Hz (e.g., model #v) prevents the model from observing any EIS magnitude saturation or the phase deviation toward zero. Under these conditions, R P is estimated as a very large resistor and the extrapolation will be different. The EIS magnitude is extrapolated from linear slope 1 and the phase will be capacitive with a value around -1Rad.
Note that α is estimated around 0.62 in the models that do not observe the low-frequency mechanism (i.e., models #iii, #iv, and #v). This is the critical α at which the temporal fitting method lost its sensitivity to R P (as was shown in figures 2(b) and in section 3.2).

EIS and the double-compartment cole model
Model #vi in figure 6 corresponds to a two-compartment Cole model. This model has sufficient degrees of freedom to provide a good fit to the EIS data in a wide frequency range. It easily reproduces the phase bump, the two linear slopes, and the quasi-saturation area between them. The two compartments of the model are expressed by {α 1 ,R P1 ,C P1 } = {0.62, 65kΩ, 2.8 μF * } and {α 2 ,R P2 ,C P2 } = {0.89, 1.5MΩ, 4.7 μF * }. Figure 7 shows the impedance spectra of these compartments superimposed with the EIS data and the total impedance spectrum of the model. The second compartment, with the higher R P , corresponds to the low frequency mechanism. The combination of α 2 and C P2 produce a linear slope and a negative phase value around -1.05 Rad at low frequencies. α 2 should be big enough to generate a large negative phase (Note: the most negative phase of a CPE may reach to / ap -2 Rad). The first compartment provides a saturated impedance (R P1 = 65kΩ) at low frequencies, then linearly pulls down the spectrum magnitude at a gradient that is mostly determined by α 1 = 0.62. The fractional order of the second compartment (α 2 = 0.89) has a lesser effect on the linear slope 1 because the impedance magnitude of the second compartment is much lower than the first compartment in mid-to-high frequencies.

Multi electrode fitting: Single versus double compartment Cole models
Inspired by the conceptual results presented in section 3.3.1. on a single electrode, we fit 4 types of Cole models to 10 electrodes. The first three models are single-compartment Cole models, and the fourth is a two-compartment model. In the first and fourth models we set f min to 0.05 Hz (the lowest measured EIS frequency), while in the second, and third models we set this parameter to 5Hz and 100 Hz, respectively (see table 2). The 5Hz cutoff frequency is large enough to disregard the low frequency mechanism, but low enough to observe the right side of the bump and the magnitude saturation. Studying this second model is beneficial as it provides a good fit to the EIS in mid and high frequencies with a lower number of parameters than the two-compartment model. Therefore, if the low frequency characteristics of an EEI is not in the scope of a study, this model is a good tradeoff between low complexity and good performance. Finally, the model with the 100Hz cutoff frequency is a good example for the models that do not have good visibility to low frequency mechanisms at the EEI. The mean and standard deviations of the fit parameters of the models are shown in table 2. Figure 8 illustrates the mean and percentiles (25%-75%) of the measured EIS data (on 10 electrodes) and the estimated EIS by the fit models. Model 1 underestimates the value of α compared to Model 2, Model 3, and the first compartment of Model 4. As explained in section 3.3.1.1., the Faradaic resistor R P in Model 1 is estimated much higher than Models 2 and 4. Model 2 underestimates the R P as it does not entirely observe the low-frequency mechanism. It partially observes the quasi-saturation part of the EIS and extrapolates from this saturation with a low R P value. Model 3 only observes the linear part of the EIS, but not the quasisaturation. The resulting extrapolation is quite different. Model 4 estimates the values of R P1 and R P2 on the order of a few MΩ and hundreds of kΩ, respectively. In addition, the second fractional order is slightly higher than the first fractional order, indicating that the linear parts have two different slopes. These factors provide much better fit for Model 4 compared to Model 1 which only has one CPE compartment. The single Cole model (i.e., Model 1) slightly deviates from the EIS at low-to-mid frequencies, whereas the double Cole model fits well to the EIS in the entire frequency range. Last but not least, in all models, little variation is observed in the estimated values of C P and R S , which stabilize around 3μF * and 500Ω respectively.

Discussions
In this study we modeled the CI EEI using the Basic and Cole models. We followed two different approaches to fit the models: A temporal approach using SPVR data, and a spectral approach using EIS data. Our experimental design was aimed at finding the factors that cause large literature discrepancies in terms of reported model parameter values.
We first investigated the temporal fitting method and its sensitivity to the double layer parameters in the Cole model (i.e., α, R P , and C P ). Investigations revealed that the region of solution is quite extended.  If the acceptance cost threshold is set to 2%, then the fractional order easily varies between 0.6 and 1.0. This is a relatively big range for α in response to such a small cost variation which is expected to be observed in literature due variabilities in measurement devices, signal filtering and processing methods, fitting methods, in vitro versus in vivo experiments, electrolyte preparations, etc. When the fractional order is set to 1.0, the temporal method estimated the values of R P and C P around 2kΩ and 65nF, respectively. This finding supports the values reported by Lim et al [8] (2-6kΩ, 10-25nF), Vanpoucke et al [9] (5-7kΩ, 5nF), Tykocinski [10] et al (2.5-5.5kΩ, 4.5-9nF), and Di Lella et al [13] (7kΩ, 6nF). However, small differences exist between the reported values in the abovementioned studies. We believe that these differences mostly stem from variations in data acquisition conditions, data processing, and model fitting methods in these studies. In addition, morphology, metal composition, roughness, and size of CI electrodes vary from one company to another. For instance, the electrode surface areas of Advance Bionics HiFocus, Cochlear Contour, and Oticon Medical EVO electrodes are around 0.2, 0.3, and 0.46mm 2 . This surface area may be used to normalize the utilized currents and to compare the current/charge densities between studies.
Although material and geometrical differences of the electrode may give varying results, we expect the presented model to be applicable to a range of electrode types.
Our temporal fitting method revealed that C P increases exponentially with the reduction of the fractional order. At low stimulation current values (<250 μA) the value of α and C P is about 0.75 and 1000nF * , respectively. This C P value is 15× greater than that of a Basic model. This result is inline with the temporally fitted α, and C P values in Mesnildrey et al [15] on CI patients which is reported to be approximately 0.64 and 481nF * respectively. Note: their electrode surface is 2.3× smaller. In addition, Mesnildrey et al estimated the value of R P to be >10 15 Ω. Our analysis showed that the temporal fitting method may result in this very high impedance value at a critical α value around 0.63. Putting these together, we may conclude that the global minimum of the cost function in Mesnildrey et al study is located at this critical α value where R P approaches infinity (the same as model #1 in figure 5).
Our time-frequency comparisons revealed that the SPVR fitting method is blind to low frequencies. Therefore, we argue that the low-value R P (i.e., a few kΩ) reported in studies like [8][9][10], and [13] should not be considered as the Faradaic resistor. These reported R P components do not have physical interpretations. The fitting cost is simply reduced by the counteracting role to the excess value of α. Precise estimation of the Faradaic resistor as well as the CPE parameters need model validity at low frequencies. The validity of a temporally fit model is 500 Hz, much higher than functioning frequencies of the Faradaic process. Even CPE characteristic is not fully observed.
Fitting Cole models using EIS data showed that more accurate results are obtained, especially at low frequencies. The fractional order was no longer estimated around 1.0 because the spectral fitting method could observe the linear slope and the negative phase in the impedance spectrum (i.e., the semi-blind area to the temporal method). As a result, the value of C P was estimated much higher (∼3.2μF * ) than those in the Basic model (∼65nF). The value of C P is estimated around 0.2μF * and 2μF * in Mesnildrey et al [15] and in Jiang et al [16], respectively. Note that in these studies, the electrode surface (2.3× smaller) and the estimated fractional order are slightly different to ours as well.
In our spectral fitting method, the fractional order was estimated to be around 0.62. It was the critical α value at which the temporal fitting method almost lost its sensitivity to R P . This value is around 0.8 in Mesnildrey et al and Jiang et al studies. We increased the lower EIS frequency band from 0.05Hz to 10Hz and then to 200Hz (i.e., the values used in those studies) but did not observe any increase in the α estimation. We hypothesize that this difference might be related to charge transfer mechanisms taking place in these studies, especially because our temporal fit showed that the value of α may change with stimulation amplitude and thus the polarization level at EEI.
Our EIS data indicates the existence of an additional charge transfer mechanism that becomes more visible in low frequencies below ∼3Hz. The co-existence of platinum and iridium materials in the electrode alloy may contribute to the complicated nature of this mechanism. Differences between the reaction potentials of these materials is expected to become more visible at low frequencies when enough time is provided for the reaction kinetics. These reactions may aid in explaining the aforementioned 'Lin. slopes' and the 'phase bump'.
Regardless of the underlying mechanism resulting in the low frequency features, it challenges the fitting of the single-compartment Cole model. Such a model does not have enough freedom of parameters to express this mechanism. In literature, several modeling configurations have been suggested for the inclusion of these effects in the EEI model. Some use a Warburg element Z W in the modeling (i.e., a CPE element with α=0.5). This element may or may not be accompanied by a parallel resistor R W . The Warburg element has been used in series with R P resistor [22,23], or as a special case of a Cole compartment (in series with R S ) [11,24,25]. We followed a similar approach presented in [26] and [24] and put two Cole compartments in series with R S to take into account the different materials (i.e., Pt and Ir) in the alloy. The Warburg element with its fixed α=0.5 could not provide the required gradient (i.e., 'Lin slope 2') and the negative phase value presented in our EIS data at low frequencies below 3 Hz, but the two compartment Cole model could provide the fit.
It should be mentioned that the validity of the twocompartment Cole model presented in this study is not guaranteed for frequencies below 0.05Hz (i.e., the lowest EIS frequency in our experiment). Working with lower frequencies may require the consideration of additional mechanisms through the use of Warburg or Cole structures in the model [24]. Other non-linearities such as those related to diffusion of reactive chemical species at low frequencies and limited thickness of the double layer [22] may need to be modeled separately. The investigation of all low-frequency mechanisms for the CI EEI is a challenging topic that needs to be addressed in future studies.
As a use case for spectral fitting with good low-frequency performance, we may refer to CI stimulation safety analysis. High charge densities or redox reactions taking place at the EEI may damage the cochlea tissue. Designing protective circuits based on a temporal fitting with underestimated values of R P and C P may not work optimally. A double (or multi) compartment Cole fit provides a better approach for the design as it can better characterize Faradaic charge transfers at the EEI.
Finally, our study suggests to clinical researchers that the parameters of electrode-tissue interface are not robustly identified using a temporal fit that is based on short stimulating pulses. If the fit model does not properly explain all frequency characteristics of the double layer, some valuable clinical effects may not be observed. In addition, a fit model that is not very sensitive to the fractional power will highly likely underestimates the values of C P and R P with a large standard deviation. Therefore, the fit parameters are not well suited to clinical and physical changes with time.
Our study also suggests that clinical researchers should observe EIS data down to very low frequencies if they are interested in finding the correct value of R P , otherwise the estimated value does not represent the Faradaic charge transfer at the electrode-tissue interface.
We believe that an in vivo study using the approach in this work would provide valuable information on the applicability of our model to the electrode-tissue interface of the cochlea. This would aid in the understanding of parameters and effects at this interface and would additionally allow for assessment of potential tissue damage [27]. Unfortunately, current CI implant technologies do not easily permit EIS data recording in vivo. Utilization of a system with sufficient sampling rate and application of the voltametric method [28,29] would enable such recordings chronically.

Conclusions
Our study conveys the following take home messages that can be considered in clinical and/or industrial studies aiming to use Cole or Basic models to evaluate a clinical condition, measure safety, or optimize a stimulation pulse or electronic circuit: (i) The temporal fitting method using SPVR data is less sensitive to variations of parameters in the Cole model than the spectral fitting method; (ii) The region of solution is very extended in the temporal fitting method (i.e., low sensitivity to model parameters). It may explain literature discrepancies and ambiguities in parameter estimations; (iii) In the temporal fitting, fixing the value of fractional order α decreases importantly the standard deviation of the estimated parameters C P and R P ; (iv) Stimulation amplitude and duration impacts the value of α; (v) Temporal fitting is blind to low frequencies so it should not be used if the lowfrequency characteristics (500Hz) of a CI EEI are in the focus of a study (e.g., long term charge accumulation across EEI in response to unbalanced stimulations); (vi) A temporally fitted Cole model using short SPVRs is only valid at mid-to-high frequencies (500Hz). This puts limitation on physical interpretations of fit parameters in clinical studies. Caution should be used when making statements about physical interpretations of a CI EEI based on the outcome of a model; (vii) The spectral fitting method extends the validity of the fit parameters to very low frequencies. However, when investigating low frequencies new mechanisms may emerge and become dominant at the CI EEI. As a result, more elements (more complexities) need to be added to the models; (viii) Enough low frequency EIS data should be observed (it is the type of investigation that determines this limit). Erroneous predictions may be obtained from a spectral-based fitted model if extrapolation of its results beyond the observed frequencies is used; (ix) Our study showed that a single-and a two-compartment Cole model may correctly reproduce the EIS data down to 5Hz and 0.05Hz, respectively.