Nonlinear effects at the electrode-tissue interface of deep brain stimulation electrodes

Objective. The electrode-tissue interface provides the critical path for charge transfer in neurostimulation therapies and exhibits well-established nonlinear properties at high applied currents or voltages. These nonlinear properties may influence the efficacy and safety of applied stimulation but are typically neglected in computational models. In this study, nonlinear behavior of the electrode-tissue interface impedance was incorporated in a computational model of deep brain stimulation (DBS) to simulate the impact on neural activation and safety considerations. Approach. Nonlinear electrode-tissue interface properties were incorporated in a finite element model of DBS electrodes in vitro and in vivo, in the rat subthalamic nucleus, using an iterative approach. The transition point from linear to nonlinear behavior was determined for voltage and current-controlled stimulation. Predicted levels of neural activation during DBS were examined and the region of linear operation of the electrode was compared with the Shannon safety limit. Main results. A clear transition of the electrode-tissue interface impedance to nonlinear behavior was observed for both current and voltage-controlled stimulation. The transition occurred at lower values of activation overpotential for simulated in vivo than in vitro conditions (91 mV and 165 mV respectively for current-controlled stimulation; 110 mV and 275 mV for voltage-controlled stimulation), corresponding to an applied current of 30 μA and 45 μA, or voltage of 330 mV at 1 kHz. The onset of nonlinearity occurred at lower values of the overpotential as frequency was increased. Incorporation of nonlinear properties resulted in activation of a higher proportion of neurons under voltage-controlled stimulation. Under current-controlled stimulation, the predicted transition to nonlinear behavior and Faradaic charge transfer at stimulation amplitudes of 30 μA, corresponds to a charge density of 2.29 μC cm−2 and charge of 1.8 nC, well-below the Shannon safety limit. Significance. The results indicate that DBS electrodes may operate within the nonlinear region at clinically relevant stimulation amplitudes. This affects the extent of neural activation under voltage-controlled stimulation and the transition to Faradaic charge transfer for both voltage- and current-controlled stimulation with important implications for targeting of neural populations and the design of safe stimulation protocols.


Introduction
Accurate representation of the electrodeelectrolyte/tissue interface is critical for the development of computational models of neural stimulation that can be used to enhance the efficacy and safety of the electrode and stimulation protocol.Computational models of deep brain stimulation (DBS) provide a valuable means with which to understand the distribution of current and voltage in the tissue surrounding the electrode and are used both clinically and in research to predict the regions of neural tissue activated during stimulation.The strength and distribution of the electric field in the vicinity of the stimulation electrodes depends on the applied current or voltage, the properties of the electrode-tissue interface, which mediates the transition of electron flow in the electrode to ion flow in the tissue, and the electrical and geometrical properties of the surrounding tissues.At the electrode-tissue interface an electrical double layer is formed, comprised of the inner Helmholtz layer with adsorbed ions and solvent molecules and a more diffuse layer of solvated ions and solvent molecules [1,2].Charge carried by electrons at the electrode is converted to charge carried by ions in the tissue either through capacitive coupling, with charging and discharging of the electrical double layer, or through Faradaic reactions involving oxidation and reduction, representing non-Faradaic and Faradaic charge transfer respectively [3].During stimulation there may also be a reduction in electrode impedance due to stimulationinduced changes in the electrode-tissue interface, including desorption of protein ions at the electrode surface [4].
At high levels of applied current or voltage, the electrode interface exhibits nonlinear behavior, with impedance decreasing due to the combined reduction in charge transfer resistance and increased double layer capacitance [2,5,6].This is accompanied by an increase in the proportion of charge transferred through Faradaic processes [7].Changes in electrode interface properties and the transition to Faradaic charge transfer, as electrode capacitance is shunted, result in a shift to the majority of charge being transferred through the resistive component of the interface, leading to reversible or irreversible chemical reactions, with the potential to cause undesirable tissue damage [8][9][10].Under voltage-controlled stimulation, these changes in the electrode-tissue interface will also affect the voltage distribution in the surrounding tissue and consequently may influence neural activation.The results of previous in vivo experimental studies in cats suggest that the current linearity limit may lie within the operating range of clinical DBS and that the electrode-tissue interface may operate within the nonlinear regime during therapeutic stimulation [7].This implies that the assumption that the DBS electrode is ideally polarizable may not always hold for clinically relevant stimulation parameters.
In the majority of computational models of DBS to-date, it is assumed that the stimulating electrode is partially or ideally polarized [11][12][13].Under this assumption, the major part of the charge is injected through the non-Faradaic component with little, or negligible, charge transferred through the Faradic component and interface properties are independent of frequency and the overpotential between the electrode and the tissue [14,15].Thus, the electrodetissue interface is typically simulated using an equivalent circuit with a capacitor or constant phase element in parallel with a charge transfer resistor [14,[16][17][18], neglecting the nonlinear and overpotentialdependent properties of the interface [19].
Based on experimental data from Richardot and McAdams [16], Cantrell and Troy incorporated nonlinear interface properties of metal microelectrodes in a finite element model using a thin layer approximation including full overpotential-dependent formulations of both resistive and capacitive interfacial components [1].Their results indicate that at low amplitudes of the applied voltage, the current density at both the electrode edge and tip is more uniform than that estimated without the electrode interface.Moreover, at higher voltage amplitudes, a reduction in impedance is observed, causing the current distribution to become nonuniform, emphasizing the need to operate stimulating electrodes at the lowest possible amplitudes [1].Similarly using a thin layer approximation, Howell et al represented nonlinear dependence of the electrode interface in a DBS model as the series combination of a Helmholtz capacitance and diffusion capacitance fitted to in vitro and in vivo experimental data [20], in parallel with a charge transfer resistance estimated using the Butler-Volmer equation.Their results indicate that during voltagecontrolled stimulation with standard DBS pulse parameters, the electrode interface can be considered as an ideal polarized electrode with a nonlinear capacitance.The charge transfer resistance, however, was assumed to be infinite which may influence the point at which the model transitions from the linear to nonlinear regime.
Nonlinearities in the behavior of the electrode interface during current-controlled stimulation have not yet been incorporated in neurostimulation models.Changes in the impedance of the electrode-tissue interface directly affect the distribution of the electric field and extent of neural activation during voltagecontrolled stimulation.While electrode impedance would be expected to have a negligible effect on the activation of neurons under current-controlled stimulation, as the voltage distribution in the surrounding tissue remains unchanged [4,18], activation overpotentials generated at the electrode interface may alter the electrochemical behavior which can have an important influence in determining the safety of the applied stimulation.
To address these questions, in this paper we implemented a nonlinear model of the electrode interface in silico representing in vitro and in vivo conditions based on rat DBS electrodes for both current and voltage-controlled stimulation.The model incorporated experimentally obtained values for the electrode interface and encapsulation tissue electrical properties [4].The unknown potentials at the electrode interface were solved using the time harmonic quasi-static formulation of Maxwell's equations with an iterative approach.The transition between linear and nonlinear behavior of the electrode interface as the applied voltage or current was increased was examined using the current-overpotential curves, or stationary polarization curves, for different frequencies.Activation of rodent subthalamic nucleus (STN) collateral axons, hypothesised to be preferentially activated during DBS [21,22], was simulated to examine the effect of the nonlinear interface on neural activation during voltage-controlled stimulation.Finally, the point of transition to the nonlinear regime during current-controlled stimulation was determined from the nonlinear polarization curves and the corresponding charge density and charge per phase were compared with the Shannon model for predicting tissue damage [23].

Methods
The relationship between the applied stimulation current or voltage and the nonlinear properties of the electrode interface of a platinum iridium electrode was examined in a model of the electrode in 0.9% saline and in a model of the rat brain.The electric field in the model of the rat brain was then coupled to a population of multicompartment neuron axon models representing branching axon collaterals within the STN to examine the effect of nonlinear electrode properties on the extent of neural activation [21,24].The electric field surrounding the DBS electrode was simulated using a 3D piecewise heterogeneous model [4].For both current and voltage-controlled stimulation, the electrode interface was incorporated into the finite element model using a thin layer approximation [16,17], with the constant phase element impedance and charge transfer resistance varying as a function of the activation overpotential in the nonlinear model.

Finite element model of DBS electrode under in vitro and in vivo conditions
A geometric model of concentric DBS electrodes was generated based on the SNEX-100 electrode (Microprobes for Life Science, Gaithersburg, USA) commonly used for preclinical rodent DBS studies.The electrode comprised an active platinum iridium contact (tip diameter 100 µm, active surface area 0.0784 mm 2 ) and a stainless-steel reference contact (diameter 330 µm, length 0.25 mm, surface area 0.34 mm 2 ).The stimulation electrode and reference contact were separated by 0.5 mm of polyimide tubing of diameter 140 µm.
To represent in vitro conditions, the concentric electrode model was placed at a depth of 1 cm in a cylindrical domain representing a plastic test tube, 2.6 cm diameter and 10 cm long, filled with 0.9% saline solution.To simulate in vivo conditions, the electrode was positioned in the STN of a heterogeneous rat brain model composed of cerebrospinal fluid, and grey and white matter.The model was created using image segmentation of the Waxholm  [28] 0.0988 164 060 White matter a [28] 0.0626 69 811 Cerebrospinal fluid [29] 2 123 Encapsulation tissue a [4] 0.05 5 × 10 5 Physiological saline [30] 1.4 143 Pt-Ir contact [11] 4.5 × 10 6 1 Steel contact [11] 9 × 10 5 1 Polyimide [11] 1 × 10 −6 1 a Values estimated at 1 kHz.
Space atlas of the Sprague Dawley rat Brain [25].
The segmented masks of the different brain tissues were converted to a geometric model using the Simpleware ScanIP software (Synopsys, USA) as described in Evers et al [4].The electrode was surrounded by 100 µm thick encapsulation tissue representing the glial scar formed during chronic electrode implantation [4].
The potential in the finite element models was estimated using the time harmonic electric-quasi static Laplace equation, where magnetic and wave propagation effects are assumed to be negligible [18,26,27], σ (ω) and ε r (ω) are electrical conductivity and relative permittivity, ω is angular frequency, ε 0 the permittivity of free space, and φ is the scalar potential.Maxwell's equation in this form considers the frequency-dependent (dispersive) conductivity and permittivity, where both conductivity and permittivity were described using the Cole-Cole equation representation of the dielectric properties of grey and white matter presented by Gabriel et al [28] and cerebrospinal fluid from De Geeter et al [29].The conductivity and relative permittivity of the encapsulation tissue were based on the values estimated from in vivo experiments in Evers et al [4].Conductivity and relative permittivity of the materials included in the models are presented in table 1.

Linear representation of electrode interface
The electrode-tissue/electrode-electrolyte interface under linear conditions was represented using the equivalent circuit model proposed by Richardt and McAdams [14,16].The circuit model consists of the parallel combination of a constant phase element, Z cpe , and resistance, R ct , with mathematical representation as given below Linear model [33] 1.42 0.85 23 Nonlinear model [1] Ae where K and β are constants denoting the magnitude of the impedance of the interface and the inhomogeneities on the electrode surface, respectively, R is the universal gas constant, F is Faraday's constant, temperature is given by T, n is the number of electrons per molecule, I o the exchange current density, i.e. the current at the equilibrium overpotential and Z EDL is total impedance of the electrical double layer interface (table 2).The constant phase element and charge transfer resistance of the electrode-tissue interface represent the non-Faradaic and Faradic charge transfer mechanisms, respectively, which are independent of frequency, activation overpotential, and the applied current or voltage amplitude in the linear model.The total admittance of the electrode interface calculated using equation ( 4) was incorporated into the in vitro and in vivo finite element models using a thin layer approximation [17,18].Due to the very thin nature of the electrical double layer, typically in the nanometer scale, mesh generation using finite element methods (FEMs) becomes computationally very expensive and can be challenging when integrating with larger structures.The thin layer approximation therefore provides a good compromise between the computational accuracy and ease of implementation [31].Following a sensitivity analysis, an exchange current density of 6.41 × 10 −4 A m −2 was assumed, based on the value used by Cantrall et al for platinum electrodes in physiological saline [17], adopted from McAdams and Jossinet [12].

Nonlinear representation of electrode interface
A similar circuit model was used to incorporate nonlinear properties of the electrode interface using parameters based on the modeling and experimental study conducted by Richardot and McAdams for platinum electrodes immersed in NaCl solution [16].
In this condition both the constant phase element Z cpe (η) and charge transfer resistance R ct (η) parameters were estimated as a function of the activation overpotential, η: where K (η) and β (η) are the overpotentialdependent form of the constants describing the magnitude of the impedance of the electrical double layer and the inhomogeneities on the electrode surface, respectively, and α a and α c are transfer coefficients.Equation ( 6) derives from the Butler-Volmer equation.Parameter values of I o = 6.41 × 10 −4 A m −2 , α a = 0.5, α c = 0.5, T = 298 K were assumed.The nonlinear constant phase element parameters were derived from the equations below as described by Cantrell et al [17], assuming that K and β are a function of the activation overpotential, η .7)  [16].A = 1.42 Ωm 2 s -β was experimentally estimated for the Pt-Ir SNEX-100 electrode in vitro [4].The activation overpotential is defined as the potential difference between the metal electrode, ϕ electrode , and the adjacent tissue or electrolyte, ϕ tissue , less the equilibrium overpotential, E 0 , It was assumed that the electrode is operating about the equilibrium overpotential [1] or that the equilibrium overpotential is equal to zero [20].Equation ( 6), the Butler-Volmer equation describes the electrochemical kinematics under the assumption that the concentrations at the electrode are not substantially different to those within the bulk tissue or physiological saline and mass transfer phenomena can thus be neglected [32].Consistent with this, postmortem examination of the brain parenchyma surrounding implanted electrodes in patients [33] and rodents [4] have not shown signs of migration of electrode material into the tissue.Similar to the linear model, the total admittance of the electrical double layer was incorporated into the thin layer approximation.A summary of the linear and nonlinear electrode interface model parameters is presented in table 2. As the overpotential is a function of the potential at the electrode interface, which is itself a function of the electrode impedance, an iterative approach was used to solve for the electric potential in the tissue surrounding the electrode as described in section 2.6.

Finite element model 2.4.1. Boundary conditions
In the finite element model, the platinum iridium contact of the microelectrode core was assigned as the active terminal and the stainless-steel contact surface was assigned as ground, figure 1(B).Neumann boundary conditions were applied to the insulating parts of the electrode and the outer surface of the rat tissue in the in vivo model, figures 1(A) and (B).A point at the dorsal roots of the spinal nerves was defined as ground, figure 1(A).For the in vitro model, Neumann boundary conditions were applied to the outer surface of the cylindrical saline domain and insulating part of the electrode.For both linear and nonlinear interface models, the electrode interface was implemented using the thin layer approximation [1] on the active contact surface for both voltage and current-controlled stimulation.The intrinsic properties and thickness of the electrical double layer were represented through the total admittance calculated from the equivalent circuit parameters in equation ( 4) which were incorporated into equation (10).
For constant current stimulation, the thin layer approximation was implemented in COMSOL using a frequency domain contact impedance boundary condition.The electrode terminal was positioned at the center of the inner core of the active metal electrode.The contact impedance boundary conditions were solved for the unknown potentials at the interface boundary with the metal, V 1 , and tissue or electrolyte, V 2 .At the interface, the normal component of the current density in domain 1 (metal), ⃗ n.J 1 , is equal to that in domain 2 (tissue), ⃗ n.J 2 : where represents the activation overpotential, η, and the admittance Y EDL is the inverse of the estimated impedance, Z EDL .For voltage control stimulation, distributed impedance boundary conditions were applied at the electrode interface in COMSOL.
The only unknown in this case is the voltage in the tissue or electrolyte, V 2 , as V 1 is the applied stimulation voltage.The normal component of the current density is again constant at the interface of adjacent domains:

Simulation of DBS waveform and neural activation
The finite element model and extracellular potential were coupled to a model of pyramidal tract axons with branching STN collaterals to examine the effect of nonlinear interface properties on activation of neurons in the region surrounding the DBS electrode.
The axon model was based on Kang and Lowery [21] incorporating membrane parameters from a mouse cortical axon model developed by Foust et al [34].The model enabled estimation of the activation threshold of the surrounding axons and collaterals incorporating electrode geometry, nonlinear electrode-tissue interface, and heterogeneity and electrical brain tissue properties.

Solver and postprocessing
The piecewise heterogeneous rat brain model with DBS electrode was meshed using Simpleware ScanIP (Synopsys), while the simpler in vitro model was meshed using COMSOL Multiphysics (COMSOL, USA).The in vivo and in vitro finite element models consisted of 5 million and 0.8 million tetrahedral elements, respectively.A quadratic interpolation function was used on each tetrahedral element to approximate the scalar potential for both models.When generating the mesh, a growth rate was assigned such that the element size increased from the smallest elements in the region surrounding the electrode and STN to a coarser mesh further from the region of interest.
Meshes with a range of resolutions were generated, from an extremely fine to a coarse resolution.The percentage change in the electric field norm and the second derivative of the electric potential over points, lines, and surfaces of the FEM model were calculated for each mesh and compared with the finest mesh.
The criterion for considering a mesh for further simulation was that the mean difference between the two must be less than 2%.The discretized finite element model was solved in COMSOL Multiphysics with a Multifrontal Massively Parallel Solver (MUMPS) direct solver using the nonlinear Newton method.Under current-controlled stimulation, with two unknowns, the potentials at the metal electrode and tissue were solved iteratively with the same solver.Application of the direct solver facilitated convergence of the solution without comprising accuracy.Initial conditions for the unknown overpotential were taken from the linear interface model to increase the speed of convergence.The stopping condition of the iterative solvers was that the relative error was less than 10 −3 .
To examine the effect of frequency on the currentand voltage-overpotential curves and the limit of linearity, simulations were first conducted at single frequencies of 130 Hz, 1 kHz, 3 kHz and 5 kHz.These frequencies were selected to represent the fundamental frequency and median frequency of the power spectrum of a typical DBS pulse of duration 60 µs-240 µs and 130 Hz stimulation frequency.Applied currents in the range −150 µA to +150 µA, were investigated for current-controlled stimulation, and voltages in the range −1 V to +1 V for voltagecontrolled stimulation.
The Fourier method was used to compare extracellular waveforms simulated using the linear and nonlinear interface models for DBS pulses with a stimulation frequency of 130 Hz, similar to those used clinically.The model was solved at the fundamental frequency of the stimulation pulse train, 130 Hz, and at 800 harmonics to span the frequency spectrum of the signal.The time domain solution was then reconstructed using the inverse Fourier transform.It should be noted that the governing equation forces all potentials in the model to vary in a sinusoidal manner [17].As a result, the solution acquired using the nonlinear impedance model at a given frequency represents the fundamental frequency only without consideration of higher order harmonics arising from nonlinear properties.In addition, the solution represents an approximation based on the weighted contribution of the constituent frequency components incorporating the nonlinear electrode properties at that frequency.The total computation time required to solve the piecewise heterogeneous rat brain model and in vitro model for a single frequency was approximately 1380 s and 30 s, respectively.Total physical memory used by the in vivo model was ∼60 GB with 2 × 2.7 GHz Intel R Xeon R processors on a Windows 2012 server.
The nonlinear region was defined as the region within which the polarization curve deviated by more than 10% from that obtained using the linear model.The root mean square difference (RMSD) between the extracellular voltage waveforms (measured at a point 1 mm from the electrode) computed using the linear and nonlinear interface models was estimated for a square pulse with pulse duration 60 µs, 120 µs, 180 µs, 300 µs and 400 µs and for an applied voltage varying from 0.1 V-5 V.
where t is the sample time, T is the total time period, y nl (t) the voltage waveform from the nonlinear interface model and y l (t) the voltage waveform from the linear interface model.The RMSD of the extracellular stimulus waveforms were not computed for currentcontrolled stimulation as the interface impedance has negligible effect on the voltage waveform under this condition.

Influence of exchange current density on polarization curves under current and voltage-controlled stimulation in the in vitro model
The relationship between applied current and overpotential was similar for both the linear model with direct solver and the iterative nonlinear model when the magnitude of the applied current lay below 45 µA for simulated in vitro current-controlled stimulation at 1 kHz, figure 2(A).When the amplitude of the applied current was increased beyond 45 µA the polarization curve obtained from the nonlinear model entered the nonlinear region with further increases in the applied current having little effect on the magnitude of the activation overpotential which remained at approximately 165 mV.A similar transition to nonlinear behavior was observed during voltage-controlled stimulation when the magnitude of the applied voltage was increased above 330 mV, with the magnitude of the activation overpotential remaining approximately constant at values beyond this, figure 2(B).
The sensitivity of the polarization curves describing the relationship between the applied current or voltage and activation overpotential to the assumed exchange current density was then examined in the in vitro simulation model for a 1 kHz source.Current-activation overpotential, I app −η, and voltage-activation overpotential, V app −η, curves were generated for exchange current densities of 1 × 10 −6 A m −2 , 6.41 × 10 −4 A m −2 , 0.01 A m −2 and 1 A m −2 .Values of the exchange current density between 1 × 10 −6 A m −2 and 0.01 A m −2 resulted in almost identical polarization curves for current and voltage-controlled stimulation.This range incorporates the typically assumed exchange current density

Influence of frequency on in vitro polarization curves
The influence of frequency on the in vitro polarization curves for the nonlinear interface model was then examined.The transition to nonlinear behavior occurred at lower values of activationoverpotential as frequency increased, at approximately 45 µA for current-controlled stimulation and 330 mV for voltage-controlled stimulation, figure 3.During current-controlled stimulation, the slope within the linear region also increased with frequency, reflecting frequency dependent changes in the impedance of the pseaudocapacitve element, figure 3(A).During voltage-controlled stimulation, a voltage divider effect is observed with the activation overpotential scaling with applied voltage within the linear region, figure 3(B).

Comparison of polarization curve for simulated in vitro and in vivo conditions
The polarization curves for the simulated in vitro and in vivo models were then compared at a frequency of 1 kHz.For both current and voltage-controlled stimulation, the electrode interface transitioned from the linear to nonlinear region at lower overpotentials in the simulated in vivo model than in the in vitro model, figure 4. For example, under current-controlled stimulation, the polarization curves entered the nonlinear region when the magnitude of the applied current exceeded approximately 30 µA with a corresponding activation overpotential of 91 mV for simulated in vivo conditions compared with 45 µA, with generated activation overpotential of 165 mV for the simulated in vitro model, figure 4(A).Under voltagecontrolled conditions, an increase in the slope of the linear region (a)-(b) was observed for the simulated in vivo conditions due to the higher impedance of the bulk tissue yielding a lower voltage across the interface than in saline for the same applied voltage, figure 4(B), similar to a voltage divider effect.An applied voltage of ±330 mV yielded a ±110 mV activation overpotential and deviated towards the nonlinear region for the simulated in vivo condition, whereas in vitro the same applied voltage amplitude yielded a much higher activation overpotential of ±275 mV at the point of transition to nonlinearity.

Influence of nonlinear electrode interface properties on the DBS voltage waveform
To examine the effect of nonlinearity of the electrode interface on the extracellular stimulation waveform during voltage-controlled DBS, the RMSD between the stimulation waveforms detected at a point inside the STN, 1 mm from the electrode, simulated using the linear and nonlinear interface models was estimated.The percentage RMSD increased exponentially as the applied voltage was increased from 0.01 V to 5 V.At the maximum voltage of 5 V the RMSD increased from 28% to 45%, 61% and 61% for pulse durations of 60 µs, 180 µs, 300 µs and 420 µs, respectively, figure 5(A).An example of an extracellular waveform simulated for linear and nonlinear interface properties for voltage-controlled stimulation with a 5 V rectangular stimulation pulse of duration 60 µs and 180 µs is shown in figure 5(B).

Influence of nonlinear electrode interface properties on the extent of neural activation
The effects of nonlinearities of the interface properties on the spatial extent of neural activation were then examined.Results are presented for voltagecontrolled stimulation only as the effect of electrode impedance on the extent of neural activation is negligible under current-controlled stimulation [4,35] (supplemental figure 1).A comparison of the electric field distribution for the nonlinear and linear model is shown in figure 6.An increase in the percentage of activated collaterals was observed when the interface entered the nonlinear region due to the associated reduction in interface impedance, figure 7. Collaterals in the STN were progressively activated with increasing pulse amplitude but at different rates for both conditions.When applying a cathodic stimulation pulse of 60 µs followed by a 120 µs charge balancing pulse of half the stimulation pulse amplitude, increasing the stimulation amplitude from 0.5 V to 0.7 V and 0.9 V resulted in activation of 0%, 14.2% and 97.3% of neuron axon collaterals for the nonlinear interface model, compared with 0%, 9.8% and 84.1% of collaterals for the linear model, (supplemental figure S2).Similarly, applying a 120 µs cathodic stimulation pulse followed by a 240 µs charge balancing pulse with amplitude 0.5 V and 0.6 V resulted in activation of 29.1% and 97.9% of collateral afferents in the nonlinear interface model, compared with only 21.7% and 81.2% of afferents in the linear model (supplemental figure S3).The difference between the two models was more pronounced for longer duration pulses, figure 7.

Influence of nonlinear electrode interface properties on the charge limits for safe stimulation
Finally, the charge and charge density at the point of transition from linear to nonlinear behavior of the interface were estimated and compared against the safety limits for electrical stimulation as predicted by the Shannon model [23].Charge and charge density were calculated for a fixed pulse duration of 60 µs with varying amplitude in the range 10 µA-1 mA.Based on the nonlinear polarization curves obtained at 1 kHz, linear operation of the interface was maintained up to a maximum pulse amplitude of 45 µA which corresponded to a maximum charge density of 3.82 µC cm −2 and charge of 3 nC when simulated using the in vitro model, figure 8. Beyond this, the model predicted that the interface operated in the nonlinear region and began to approach the FDA approved limit for safe charge transfer of 30 µC cm −2 .The 30 µC cm −2 limit corresponded to an amplitude of 360 µA and charge of 23.4 nC for the pulse durations examined.As observed in the polarization curves, the transition to nonlinear behavior occurred at slightly lower values of 2.29 µC cm −2 and 1.8 nC in the in vivo model, figure 8(B).

Discussion
Nonlinear behavior of the electrode interface impedance at high values of current or voltage are wellestablished in experimental studies.However, computational models rarely if ever incorporate these nonlinearities.Incorporating the nonlinear interface is particularly complex for current-controlled stimulation where both the potential generated at the electrode and the potential at the tissue are unknown.In this study, nonlinear interface properties of DBS electrodes were incorporated in a finite element model using a thin layer approximation and iterative solver for both current and voltage-controlled stimulation under in vitro and in vivo conditions.The model captured both the established nonlinearity of the interface and frequency dependency of the polarization curves.Inclusion of nonlinear properties affected collateral activation during voltage-controlled stimulation, with the difference between linear and nonlinear models increasing with pulse duration.Although the extent of neural activation was not affected by interface impedance under current-controlled stimulation, the transition to Faradaic charge transfer at higher applied currents with potential reversible and irreversible electrochemical reactions at the electrode surface may result in electrode or tissue damage.The model simulations predict that this transition to Faradaic charge transfer occurs at current densities well-below the Shannon safety criteria and within the therapeutic range of clinical DBS.

Nonlinear model and sensitivity to exchange current density
The simulated polarization curves from the linear interface model overlapped with the nonlinear model when the magnitude of the applied current or voltage was less than approximately 45 µA or 330 mV, respectively, for a 1 kHz source, figures 2(A) and (B).With the linear model, the electrode is effectively capacitively coupled to the tissue with current transfer primarily occurring through charging and discharging of the interface with negligible Faradic current [20,36].However, when the magnitude of the applied voltage or current increased beyond those limits, the charge transfer mechanism depended on the activation overpotential generated by the applied current or voltage.A sensitivity analysis of the influence of the exchange current density, confirmed that the polarization curves for both current and voltage-controlled stimulation were relatively insensitive to the value of the exchange current density within the experimentally reported range for platinum and Pt-Ir [14,16,17,[36][37][38], figures 2(C) and (D).This is consistent with Richardot et al who showed that variation in the exchange current density was not the main source of nonlinearity at the electrode but acts as a multiplying coefficient rather than changing the current potential characteristics [16].Examining simulated  metal implants in the vicinity of electrical stimulation, Mercadal et al observed that the bounds on the transition region where the linear relationship between electric field and current density through the electrode breaks down, was also insensitive to exchange current density [38].

Frequency-dependency of polarization curves and effect of volume conductor electrical properties
The influence of frequency on the polarization curves was apparent for both current and voltage-controlled stimulation, figure 3.Under voltage-controlled stimulation, the influence of frequency on the linear region of the polarization curve was negligible.When the applied voltage was sufficiently low (<200 mV), similar activation overpotentials were observed in the range examined between 130 Hz and 5 kHz.However, when the magnitude of the applied voltage increased beyond approximately 330 mV, the activation overpotential decreased with increasing frequency, and nonlinear interface behavior was observed at lower activation overpotentials as frequency increased, figure 3(B).This corresponded to higher values of the interface current at higher frequencies, driven by the reduction in the impedance of the constant phase element impedance, and overall interface impedance, with increasing frequency [2,16,39].The linear region of the polarization curves under current-controlled stimulation exhibited a frequency-dependency due to the influence of frequency on the constant phase element impedance and through this on the activation overpotential, figure 3(A).The point of transition to nonlinear behavior occurred at similar values of applied current within the frequency range examined.
The polarization curves differed also between simulated in vitro and in vivo conditions for both current and voltage-controlled stimulation, figure 4.Under controlled-current stimulation in vivo, the interface entered the nonlinear region at lower values of applied current and activation overpotential when compared with the simulated in vitro model, with the linear region of the polarization curves identical up to the transition point of 28 µA at 1 kHz.Under voltage-controlled stimulation, the electrode interface also entered the nonlinear region at lower values of the activation overpotential for simulated in vivo than in vitro conditions.However, this occurred at similar values of applied voltage and the slope of the linear region of the polarization curve in vivo was higher than that in vitro, due to the electrical properties of the electrolyte and tissue around the DBS electrode.Under simulated in vitro conditions, the electrode was surrounded by physiological saline which has a lower impedance than the encapsulation tissue and surrounding grey matter in vivo [40].
In the linear region, the voltage divider effect across the electrode and bulk tissue causes the activation overpotential for a given applied voltage to be lower for the in vivo than in vitro model, driving the interface into the nonlinear regime at lower voltages.This has important implications for the design of stimulating electrodes as most studies of the electrode electrochemical behavior have been conducted under in vitro conditions using cyclic voltammetry.Wei and Grill similarly observed a higher charge transfer ratio in vitro than in vivo, suggesting that the non-Faradaic charge capacity of clinical DBS electrodes is overestimated in vitro [7].Our results support this observation and identify the influence of the impedance of the surrounding volume conductor on R ct and Z cpe , through the overpotential, as a contributory factor.Differences in chemical reactions at the electrode in vivo and in vitro may also further affect the charge transfer.

Extent of collateral activation under voltage and current-controlled stimulation
When comparing the extent of activation of axon collaterals estimated using the linear and nonlinear interface models, the linear interface model underestimated the fraction of collateral branches activated at higher amplitudes for voltage-controlled stimulation.The RMSD between stimulation waveforms simulated with the linear and nonlinear interface models increased with applied voltage, figure 5(A).Howell et al [20] estimated that the RMSD of voltage waveforms estimated using linear and nonlinear interface models was less than 0.4% for a 5 V pulse of 210 µs duration [20].A similar 5 V stimulation pulse of duration 180 µs, resulted in a RMSD of approximately 45% using the model presented here, figure 5(B).The observed differences between the two models may be attributed to the presence of the nonlinear Faradaic resistance which accounted for the majority of the charge transfer in the nonlinear regime in the present model, but was not included in [20], and, to a lesser extent, on the presence of encapsulation tissue surrounding the DBS electrode in the model presented here which tends to narrow the linear operating regime (see supplemental figure S5).
Under current-controlled stimulation, nonlinearity of the interface had a negligible effect on the extent of collateral activation (supplemental figure S1).It is well established that current-controlled stimulation maintains high efficacy in stimulating neurons despite the presence of variations in impedance caused by dynamic changes at the electrode-tissue interface [41].However, the safety and electrochemical stability of the electrode operating under currentcontrolled stimulation are dependent on the activation overpotential.Limiting current injection to non-Faradaic charge transfer is, therefore, a high priority in maintaining electrochemical stability [36].For currents below 45 µA the interface remained in the linear region of the polarization curve under simulated in vitro conditions, reducing to 30 µA for simulated in vivo conditions.This corresponds to a charge density of 3.82 µC cm −2 and charge of 3 nC for a typical DBS pulse with duration 60 µs in vitro conditions and 2.29 µC cm −2 and 1.8 nC, respectively, in vivo.

Relation between region of nonlinearity and Shannon safety criteria
The results predict that for rodent DBS electrodes, the interface impedance enters the nonlinear region at currents and voltages below the Shannon limit for tissue damage of 4 nC/phase for microelectrodes and 30 µC cm −2 for macroelectrodes [23], figure 8.The transition point predicted in the presented model also lies well-within the range of stimulation parameters used in rodent models of DBS [4,42].These represent estimated values under the conditions simulated.The exact current or voltage at which the limit of nonlinearity occurs for a given electrode will depend on the properties of the surrounding tissues and the electrode interface, and the frequency of the applied input.In addition, the spatial distribution of the current density across the surface of the electrode will vary with electrode size and shape.Nevertheless, the ability to model the transition to Faradaic charge transfer provides a means to incorporate variations in stimulation frequency, pulse duration and amplitude, and electrode and surrounding tissue properties, that are not captured in the Shannon model nor included in current DBS models.In addition, it provides an estimate of the point above which the consequences of Faradaic charge transfer, such as pH changes and redox reactions, are likely to occur.In this way DBS models incorporating non-linear effects can be used to optimize stimulation paradigms to ensure safety and stability of the electrode-tissue interface for long term stimulation.Faradaic reactions such as reduction, which involves the addition of an electron and oxidation involving the removal of an electron at the electrode, are either reversible, when formed products remain bound to the electrode and can be recovered, or irreversible, when formed products diffuse away from the electrode and cannot be recovered.Within the current model, it is not possible to distinguish between reversible and irreversible Faradaic reactions.Incorporation of mass transport phenomena to capture this would enable effects such as alterations in pH and corrosion onset to be estimated.
Though developed for rodent models of DBS, the methods presented can be implemented in patient DBS models.If the cutoff for the nonlinear region is extrapolated for the larger DBS electrodes used in humans, a similar picture emerges.Scaling to human DBS electrodes with a typical surface area of 0.06 cm 2 and 60 µs pulse duration, the point at which the  [43].The estimated limit of nonlinearity obtained from the model is higher than the RMS amplitude of 0.02 mA estimated experimentally for the current linearity limit of a DBS electrode of a surface area of 0.058 cm 2 by Wei and Grill [7].However, in that study the limit was estimated for frequencies less than 100 Hz and a linearity limit was not evident at frequencies above 100 Hz in the range 0.01 mA-0.2 mA examined.
The results support the use of DBS stimulation protocols that that would enable lower stimulation amplitudes to avoid operating within the nonlinear regime, and that ensure charge recovery to avoid charge accumulation in the case of Faradaic reactions.Charge balancing of the stimulation waveform is commonly achieved through the use of biphasic stimulation waveforms, where a second charge balancing pulse is designed to deliver an equal amount of charge, of opposite polarity, to the leading stimulating waveform [36].A number of other strategies have also been proposed to avoid electrode or tissue damage, including active discharging of the interface through a capacitive circuit.In addition to ensuring charge recovery, the magnitude of the charge balanced waveform should be limited to avoid electrode corrosion due to irreversible Faradaic reactions.Longer pulse durations could potentially be used to reduce stimulation amplitude, whilst maintaining equivalent levels of neural activation.Aside from electrochemical reactions, it should be noted that non-Faradaic reactions can also cause tissue damage due to repeated depolarization or hyperpolarization of surrounding neurons [44].

Assumptions of the nonlinear electrode interface model
Although the nonlinearity of the electrode-tissue interface impedance is well-established, relatively few studies have compared linear and nonlinear interface models, particularly in the area of DBS.Cantrell et al [1] compared the effect of linear and nonlinear models of the electrode-electrolyte interface on the current density distribution, observing increasing nonuniformity of the current distribution at higher driving potentials as the profile approximated that generated when the interface was neglected.Howell et al [20], also compared linear and nonlinear models.In contrast with the present study in which the observed nonlinearity is primarily due to the Faradaic charge transfer resistance [14], the charge transfer resistance was assumed to be negligible and the electrical double layer was approximated as an ideal polarized electrode represented by a nonlinear capacitance.This results in the maintenance of a relatively high impedance across the electrode as the applied current or voltage is increased (see supplemental figure S5).When a nonlinear charge transfer resistance is included, a transition to nonlinearity is observed similar to that observed here, though the point of transition occurs at different driving currents.
The model presented is a simplification of the biophysical processes involved and has limitations which should be considered when interpreting the results.The electrochemical interface model does not consider the overpotentials induced by changes in the concentration of ions in the surrounding medium.Mass transport phenomena are thus not incorporated, and the model does not describe corrosion of the electrode and its corresponding diffusion overpotential.The transfer coefficients for anodic and cathodic and equilibrium overpotential were assumed to be 0.5 and zero, which will influence the polarization curves, though for the fixed geometric dimension of the electrode, the interface model has been shown to be relatively insensitive to the transfer coefficients unless corrosion kinetics are involved [38].While the model predicts the point of region of transition between Faradic and non-Faradic charge transfer on the polarization curves it does not distinguish between reversible and irreversible Faradaic charge transfer.The nonlinear interface parameters were assumed to be the same for in vitro and in vivo conditions.Differences between the two may also affect charge injection mechanisms.The activating overpotential is incorporated in the Butler-Volmer equation (equation ( 6)) which does not capture the nonlinearity caused by harmonics of the signal on the polarization curve.The parameters of the pseudocapacitance and Faradaic resistance were assumed to be independent of frequency.An explicit dependency of these parameters on frequency would alter the observed relationship between the limit of nonlinearity and frequency, which arose here through the relationship between pseudocapacitive impedance and frequency, and the resulting effect on the overpotential.In simulating the square wave pulses, the Fourier approach assumed the response at each frequency to be independent of other frequencies which may shift the nonlinear response to higher amplitude currents or voltages.Finally, as the real electrical double layer has a volume with electrical properties that are usually non-uniform, the mathematical representation of the thin layer should be considered as an approximate representation.

Conclusion
In this study, nonlinear properties of the electrode interface were incorporated in an finite element model of rat DBS electrodes for both current and voltage-controlled stimulation.The model predicts a transition to nonlinear behavior of the interface impedance for applied voltages or currents, at 1 kHz, above approximately 330 mV or 45 µA, respectively.For a 60 µs stimulation pulse, the limit of linearity was estimated as 2.29 µC cm −2 .The results indicate that the threshold for nonlinearity varies with experimental condition and, for currentcontrolled stimulation, with the properties of the surrounding tissue.During voltage-controlled stimulation, nonlinearity of interface had a substantial effect on the extracellular stimulation waveform and thresholds for neural activation which could be a confounding factor during DBS programming.While the extracellular potential is insensitive to changes in impedance during current-controlled stimulation, the threshold of nonlinearity corresponds to a transition to Faradaic charge transfer which, in the case of irreversible reactions, may result in electrode or tissue damage.When extrapolated to larger human DBS electrodes, the results suggest that the transition to nonlinear behavior and Faradaic charge transfer may occur within the range of clinically-relevant stimulation parameters.Inclusion of nonlinear electrode properties in DBS models provides a more complete picture of electrode-tissue interface stability during long-term stimulation and can help with optimizing electrode design for high electrochemical stability.

Data availability statement
The data cannot be made publicly available upon publication because the cost of preparing, depositing and hosting the data would be prohibitive within the terms of this research project.The data that support the findings of this study are available upon reasonable request from the authors.

Figure 1 .
Figure 1.Representation of the thin layer approximation at the electrode and boundary conditions for bipolar stimulation for (A) in vivo and (B) in vitro finite element models of the rat DBS electrode.

Figure 2 .
Figure 2. (A) Applied current-activation overpotential curve and (B) applied voltage-activation overpotential curve for simulated in vitro current-controlled stimulation and voltage-controlled stimulation, respectively.Simulations were conducted for the SNEX-100 microelectrode with an applied stimulation current or voltage at 1 kHz.Data are presented for linear (orange) and nonlinear (blue) interface properties, the transition points from the linear region to the nonlinear region (nonlinear boundary) are indicated.(C) & (D) Sensitivity of applied current-overpotential and applied voltage-overpotential curves to exchange current density, Io, in the in vitro model.The regions a-b and a ′ -b ′ represent the linearized kinematic regions under current-controlled stimulation (C) and voltage-controlled stimulation (D).

Figure 3 .
Figure 3. Current-overpotential and voltage-overpotential curves simulated for in vitro conditions at 130 Hz (yellow), 1 kHz (blue), 3 kHz (red) and 5 kHz (green).The upper and lower nonlinear boundaries are provided for each frequency (for 130 Hz in yellow, 1 kHz in blue, 3 kHz in red and 5 kHz in green).Data are presented for (A) current controlled stimulation and (B) voltage-controlled stimulation.

Figure 4 .
Figure 4. Comparison of applied current-activation overpotential and applied voltage-activation overpotential curves estimated for simulated in vivo (blue) and in vitro (orange) conditions at a frequency of 1 kHz for (A) current controlled stimulation and (B) voltage-controlled stimulation.

Figure 5 .
Figure 5. (A)Percentage root mean square difference (RMSD) between extracellular DBS stimulation waveforms simulated using linear and nonlinear interface parameters at a point inside the STN, located 0.5 mm from the electrode.The RMSD is expressed as a percentage, normalized with respect the waveform simulated using linear interface properties.(B) Voltage waveform measured 0.5 mm from the electrode simulated with the linear and nonlinear electrode interface for a square wave stimulation pulse of duration 60 µs and 180 µs and amplitude of 5 V.

Figure 6 .
Figure 6.Simulated in vivo electrical field distribution.Comparison of nonlinear (left) and linear (right) model electric field distribution around the stimulation electrode in the rat STN at 0.5 V, 2 V and 4 V voltage-controlled stimulation at 1 kHz.

Figure 7 .
Figure 7. Simulated axon collateral activation under voltage-controlled stimulation.(A) Comparison of stimulated (green) and unstimulated (red) STN collaterals in the rat for different pulse amplitudes of 0.7 and 0.9 V simulated with nonlinear electrode-tissue interface properties (left) and linear electrode-tissue interface properties (right).(A) Data were simulated for a biphasic pulse with cathodic stimulation pulse of duration 60 µs followed by a 120 µs anodic pulse.(B) Data were simulated for a biphasic pulse with cathodic stimulation pulse of duration 120 µs followed by a 240 µs anodic pulse.

Figure 8 .
Figure 8. Region of charge and charge density where linear and nonlinear behavior is observed in the model presented alongside the Shannon model for safe levels of stimulation derived for a microelectrode with surface area 7.85 × 10 −4 cm 2 with k = 1,1.5 and 1.85, under current-controlled stimulation.The charge density per phase and charge per phase were calculated for a square pulse width of 60 µs with varying current amplitude in the range of 10 µA to 1000 µA in 10 µA intervals.Solid green triangles correspond to stimulation amplitudes within the linear regime of the electrode interface model and solid red circles represent stimulation within the nonlinear regime.The green dotted line represents the maximum charge density of the electrode in the linear region for nonlinear interface and the red dotted line represents the corresponding charge per phase.(A) in vitro nonlinear model and (B) in vivo nonlinear interface model.

Table 1 .
Conductivity and relative permittivity values of the materials included in the in vivo and in vitro DBS models at 1 kHz.

Table 2 .
Parameter values for linear and nonlinear electrode interface models.

Table 3 .
Stimulation current for clinical DBS electrodes at a charge density equivalent to that at the limit of linearity predicted for the simulated SNEX-100 preclinical electrodes for a stimulation pulse of 60 µs pulse duration.