Journal of Neural Engineering

PAPER • OPEN ACCESS You may also like

- Delayed administration of interleukin-4
Choosing the right electrode representation for coacervate alleviates the neurotoxic
phenotype of astrocytes and promotes
modeling real bioelectronic interfaces: a functional recovery after a contusion spinal
cord injury
Manoj K Gottipati, Anthony R D’Amato,
comprehensive guide Jayant Saksena et al.

- Connectivity study on resting-state EEG

To cite this article: Aleksandar Opančar et al 2024 J. Neural Eng. 21 046049 between motor imagery BCI-literate and
BCI-illiterate groups
Hanjin Park and Sung Chan Jun

- A preliminary study exploring the effects of

transcutaneous spinal cord stimulation on
View the article online for updates and enhancements. spinal excitability and phantom limb pain in
people with a transtibial amputation
Ashley N Dalrymple, Lee E Fisher and
Douglas J Weber

This content was downloaded from IP address 185.80.143.159 on 08/10/2024 at 19:34

 J. Neural Eng. 21 (2024) 046049 https://doi.org/10.1088/1741-2552/ad6a8b

Journal of Neural Engineering

PAPER

Choosing the right electrode representation for modeling real

OPEN ACCESS
bioelectronic interfaces: a comprehensive guide
RECEIVED
28 February 2024 Aleksandar Opančar1,2, Eric Daniel Głowacki2,∗ and Vedran Ðerek1,∗
REVISED 1
12 July 2024 Department of Physics, Faculty of Science, University of Zagreb, Bijenička c. 32, 10000 Zagreb, Croatia
2
Bioelectronics Materials and Devices Laboratory, Central European Institute of Technology, Brno University of Technology, Purkyňova
ACCEPTED FOR PUBLICATION
123, 61200 Brno, Czech Republic
1 August 2024 ∗
Authors to whom any correspondence should be addressed.
PUBLISHED
14 August 2024 E-mail: glowacki@vutbr.cz and vdjerek@phy.hr

Keywords: bioelectronics, electrodes, interface, constant phase element, simulation

Original content from
this work may be used Supplementary material for this article is available online
under the terms of the
Creative Commons
Attribution 4.0 licence.
Any further distribution Abstract
of this work must
maintain attribution to
Objective. Producing realistic numerical models of neurostimulation electrodes in contact with the
the author(s) and the title electrolyte and tissue, for use in time-domain finite element method simulations while
of the work, journal
citation and DOI. maintaining a reasonable computational burden remains a challenge. We aim to provide a
straightforward experimental-theoretical hybrid approach for common electrode materials (Ti,
TiN, ITO, Au, Pt, IrOx) that are relevant to the research field of bioelectronics, along with all the
information necessary to replicate our approach in arbitrary geometry for real-life experimental
applications. Approach. We used electrochemical impedance spectroscopy (EIS) to extract the
electrode parameters in the AC regime under different DC biases. The pulsed electrode response
was obtained by fast amperometry (FA) to optimize and verify the previously obtained electrode
parameters in a COMSOL Multiphysics model. For optimization of the electrode parameters a
constant phase element (CPE) needed to be implemented in time-domain. Main results. We find
that the parameters obtained by EIS can be used to accurately simulate pulsed response only close
to the electrode open circuit potential, while at other potentials we give corrections to the obtained
parameters, based on FA measurements. We also find that for many electrodes (Au, TiN, Pt, and
IrOx), it is important to implement a distributed CPE rather than an ideal capacitor for estimating
the electrode double-layer capacitance. We outline and provide examples for the novel
time-domain implementation of the CPE for finite element method simulations in COMSOL
Multiphysics. Significance. An overview of electrode parameters for some common electrode
materials can be a valuable and useful tool in numerical bioelectronics models. A provided FEM
implementation model can be readily adapted to arbitrary electrode geometries and used for
various applications. Finally, the presented methodology for parametrization of electrode materials
can be used for any materials of interest which were not covered by this work.

1. Introduction [1–7]. Stimulators can be external or implanted,

wired or wireless, externally, or battery-powered [8–
The use of electronic devices in living systems for 12], whereas stimulation protocols can be open- or
both electrical stimulation and recording requires closed-loop, continuous or intermittent, depending
a wide variety of bioelectronic methods and tech- on the needs of the application. However, all mod-
niques, including single-cell patch-clamp for study- alities of electrical stimulation or electrical signal
ing cell membrane properties, microelectrode arrays recording require a pair of electrodes as interfaces
for cell culture or tissue stimulation and recording, all with the stimulation/recording target, where the cur-
the way to the central nervous system and peripheral rent enters or leaves the region. The efficient and
nerve implants for therapeutic and research purposes safe application of electric current to living cells or

© 2024 The Author(s). Published by IOP Publishing Ltd

J. Neural Eng. 21 (2024) 046049 A Opančar et al

tissues for purposes of either research or medical coupled to the stimulation and return electrodes,
treatment is governed by the properties of the stim- via the model of the environment. Several physical
ulation electrodes [13, 14]. Charge transport from and electrochemical processes may occur at the elec-
the stimulation device to the electrode is electronic trode/electrolyte interface, such as formation of capa-
in nature, whereas in tissue, it is purely ionic [15]. citive double layers or Faradaic reduction or oxid-
The electrode/tissue interface is where the transduc- ation processes, which may be further enhanced or
tion of electronic to ionic transport occurs, and the limited by the transport of redox active species in
physical properties of the electrodes, together with the vicinity of the electrode in the electrolyte. Ab ini-
the stimulation protocol, determine if the electrical tio models considering all the species in the electro-
coupling will be capacitive, Faradaic, or of a mixed lyte, as well as the electrode material properties at
nature [16, 17]. A Faradaic electrode will facilitate the atomistic scale are neither feasible or reasonable
the charge transfer by a redox reaction, either redu- with the current computational abilities, at least not
cing or oxidizing chemical species in the surround- for large spatial- and temporal-scale models. To effi-
ing environment or the electrode itself, resulting in ciently numerically simulate the electrode response to
the release of ions and free radicals, or in corrosion of a voltage stimulus in an electrolyte solution, a model
the electrode. Typically, biphasic current pulsing pro- with a minimal number of parameters is desired,
tocols are applied to maintain charge balance, how- considering its geometry, material, and electrochem-
ever, Faradaic charge transfer reactions may not be ical properties, especially in spatially extensive three-
reversible and can cause net changes in the chemical dimensional models where each spatial mesh element
environment surrounding the electrode [18]. Both has multiple degrees of freedom associated with it.
processes can result in failure in stimulation and/or This electrode can then be coupled with other subsys-
recording, as well as unwanted interactions of the tems within a comprehensive numerical simulation.
living system with species released by redox reac- We will show how said electrode properties can
tions. Therefore, when designing stimulation devices be extracted from simple electrochemical measure-
and stimulation protocols, due attention should be ments, such as electrochemical impedance spectro-
paid to the electrodes, their geometry, and their elec- scopy (EIS) with applied bias and fast amperometry
trochemical properties. Ultimately, only experiments (FA) under constant potential. Using the parameters
will conclusively determine the safety and efficacy extracted and estimated from those measurements,
of the stimulation device, protocol, and electrodes the electrode/electrolyte can be faithfully implemen-
[19, 20]. In vitro and in vivo experiments carry a ted in a three-dimensional time-dependent finite ele-
substantial cost of time and resources, as well as an ment method multiphysics model, which can be used
ethical burden, especially when iterations in design for arbitrarily shaped electrical stimulation or record-
need to be statistically verified in large populations of ing. In particular, it is vital that the model is time
laboratory model [21]. With rapid progress in com- domain based opposed to frequency domain to make
putational capabilities, the concept of digital twins it in principle compatible to a Hodgkin–Huxley type
has been widely exploited to allow purely numer- model which can be used to calculate the neural
ical experiments—simulations of a biological system response. We will compare commonly used paramet-
in a reciprocal interaction with the environment—to rizations of the electrode interface materials which
supplement and reduce the necessity for real exper- are relevant to the field of bioelectronics research
iments, thus speeding up the research and reducing (gold, platinum, titanium, indium tin oxide, titanium
its financial and ethical burden [22–24]. A compre- nitride, and iridium oxide). The electrochemical
hensive digital twin simulation should include both properties of the material will be modelled as an
the living and electronic sides of the system with a equivalent Voigt RC circuit with voltage-dependent
faithfully represented cell or tissue model, as well as double layer capacitance and Faradaic resistance, and
the full model of the stimulation device. To fully com- as the Randles circuit with voltage-dependent charge
ply with the requirements, all the digital twin model transfer resistance and constant phase element (CPE)
parameters should be dynamically updated by the real parameters distributed along the electrode surface.
measurements. Finally, all components of the simula- We will demonstrate the implementation of both
tion must be individually experimentally verified, and Voigt and Randles distributed equivalent circuits in
the applicable range of the parameters established. the time-domain within the COMSOL Multiphysics
To fully describe electrical stimulation in bio- environment. The parameters for both equivalent cir-
electronic systems, it is necessary to model the cur- cuits will be extracted from the EIS measurements,
rent injection from a bioelectronic device through the and we will show that our model faithfully reproduces
electrode, and into the electrolyte, where the current the measured EIS data in the frequency domain. We
faces three-dimensional objects which may or may will show that the parameters extracted from EIS can-
not be electrically active as well. Finally, the current not be directly applied to pulsed stimulation, due to
is extracted through the return electrode. Provided the inability of the EIS measurements to capture tran-
that a faithful device model exists, it will be directly sient phenomena that can be present in low duty cycle

2
J. Neural Eng. 21 (2024) 046049 A Opančar et al

pulsed stimulation, due to transient changes in elec- characterize its transport properties is advantageous,
trode properties or rapid depletion of redox-reactive because the low number of parameters required for
species in the vicinity of the electrode when rap- the model can be efficiently and rapidly solved using
idly brought out of the equilibrium, especially when the state-of-the-art hardware, even with a massive
measuring the EIS under bias voltage. We will demon- number of degrees of freedom used in a simulation.
strate and verify a method for self-consistent correc- There are many commonly used equivalent circuits,
tion of the EIS parameters by nonlinear fitting of but two of the simplest and most commonly used are
the equivalent circuit parameters to the fast ampero- the Voigt circuit and the simplified Randles circuit
metry measurements (FA), starting from the para- with a CPE (figure 2) [28, 29].
meters obtained by the EIS. Finally, we will present The Voigt circuit (figure 2(a)) has the advant-
the validation of our model by extrapolation to differ- age of having a very clear interpretation of each ele-
ent geometries, changed electrolyte conductivity and ment, containing an electrolyte resistance in series,
stimulation protocols. with the double layer capacitance and charge trans-
A number of recent works have tackled the fer resistance in parallel. In addition, its numerical
numerical modelling of electrodes in contact with implementation is straightforward both in the time
excitable tissues, each of which focused on a unique and frequency domains, since all elements can be rep-
aspect of the problem, such as the electrode shape resented by linear differential equations. On the other
optimization [25, 26] or the studies of the effective hand, the Randles circuit contains a CPE [30].
electrical parameters at the stimulation target [27], All equivalent circuit elements are well defined in
while some have attempted to take into account also the frequency domain by their complex impedances,
the electrochemical properties of the electrode [24] ZR = R, ZC = jω1C and ZCPE = Q(jω) 1
α , with α ranging

in the regime where no Faradaic reactions occur. from 0–1 for the capacitive CPE [31, 32]. The res-
Our attempt, however, includes a comprehensive ulting impedance can be easily evaluated in the fre-
FEM approach to the electrode impedance, using quency domain, but can pose difficulties when imple-
both distributed linear equivalent circuits and distrib- menting time-domain models.
uted CPE equivalent circuits, with all the paramet- It can be seen from the complex impedances that
ers being extracted from the measurements and the CPE falls somewhere in between the resistor and the
model further optimized to reproduce the validation capacitor depending on the exponent α, where α
stimulation pulses. This approach has enabled us to = 0 reproduces the resistor and α = 1 the capacitor
identify both the materials and the electrical poten- with the parameter Q having the role of conductance
tial windows for which a simplified linear simulation and capacitance respectively. CPE with exponent α
approach may be valid. Furthermore, due to the used between those limiting values has a rather complex
process of model validation by using the stimulation physical interpretation, but it models a behaviour
pulses, our model has inherent ability to be used as a very commonly encountered in electrochemistry.
digital twin, by feeding the measured data dynamic- The CPE behaviour of the electrode can either
ally back to the model for validation and parameter be an intrinsic microscopic property of the material,
optimization [22, 23]. or a global behaviour originating from distribution
To our knowledge, this is the first report of of reactivity due to interface inhomogeneities such
a non-approximate implementation of distributed as surface disorder and roughness, electrode poros-
electrochemical equivalent circuits with constant- ity, specific ion adsorption, electrode geometry and
phase elements in the time domain within the normal-to-surface distributions of properties in films
COMSOL Multiphysics FEM environment, as well and coatings [33]. In all the mentioned cases, it res-
as a unique comprehensive overview of electrochem- ults from a 2D variation of properties along the sur-
ical properties for commonly used electrode mater- face of the electrode or 3D variation of properties
ials. Ultimately, the electrode model presented can also in the direction normal to the electrode surface
be used to accurately represent an electrode/electro- [34]. Putting aside the microscopic properties of the
lyte interface used for stimulation or measurement in electrode which can lead to a CPE behaviour, even
time domain, within a comprehensive finite element an ideally capacitive electrode will have an apparent
method simulation. CPE behaviour above a certain limiting frequency flim
due to the finite size of the electrode and the edge
2. Method effects, leading to inhomogeneous electric field dis-
tribution across the electrode surface [35, 36]. Here,
2.1. Representation of electrochemical properties the apparent CPE behaviour refers to the fact that
of the electrodes by equivalent circuit models CPE exponent α is a function of frequency which is
Electrodes in numerical simulations involving elec- not the case for a normal CPE element. For the disc
trochemistry are typically represented by an equival- electrode of radius r the expression for flim is given
ent electrical circuit. A concept of estimating elec- by flim = 2πκCr , where κ is the electrolyte conductivity
trical equivalent circuit parameters on an electrode to and C is the specific capacitance of the electrode [33].

3
J. Neural Eng. 21 (2024) 046049 A Opančar et al

For an electrode with a radius of 1 mm, the capacit- solution. A 10× diluted version of the same electro-
ance C can span a range of two orders of magnitude lyte, 0.1x PBS, was prepared by diluting the 1xPBS in
for commonly used materials, from approximately 10 1:9 ratio with DI water. An Ag/AgCl wire was used as
µF cm−2 for indium tin oxide (ITO) to 1000 µF cm−2 a reference electrode (REF) and placed 5 mm from
for iridium (IV) oxide (IrOx). That gives the range the WE. The electrodes were connected to the cor-
of flim from 2.4 kHz to 24 Hz, well within the typical responding ports by alligator clips to a PalmSens 4
frequency range used in bioelectronics, which makes or Ivium PocketStat2 potentiostat, or alternatively to
the effects of flim crucial for consideration of these a combination of a potentiostat and a digital oscillo-
applications. scope (Pico Technology, Picoscope 4424 A).
In our work, we will use both Voigt and Randles
equivalent circuits to represent electrochemistry on 2.4. EIS
the electrode surface, to evaluate the optimal condi- EIS was measured in a 3-electrode setup, using the
tions in which they may be used. Further, we will con- phosphate-buffered saline (PBS) as an aqueous elec-
sider all of the equivalent circuit parameters to be in trolyte and a platinum CE. Measurements were con-
principle voltage-dependent. ducted with PalmSense4 potentiostat at a DC voltage
bias VP , ranging from −900 + 800 mV in steps of
100 mV vs. Ag/AgCl reference wire, in a frequency
2.2. Electrode preparation
range of 100 mHz–100 kHz with 8.5 frequencies per
Microscope slides (1 × 1) inch2 were cleaned by
decade (52 frequencies in total per DC bias.) The AC
successive ultrasonication in acetone, isopropanol,
amplitude was set to 10 mV for all measurements.
Micro-90 detergent, and DI water, and finally treated
with oxygen plasma (Diener NANO Plasma Cleaner).
2.5. FA
All samples were then sputter coated with a 100 nm
FA was measured in the same electrochemical cell and
layer of Ti using a Kaufman ion-beam source (IBS).
under the same conditions as the EIS measurements.
This Ti acts as a common layer below all studied
Current response to the constant voltage pulses
samples with the exception of ITO, as it has excel-
applied to the WE, ranging from −900 + 800 mV
lent adhesion on glass and is a suitable underlayer for
in steps of 100 mV vs. Ag/AgCl was measured We
all the studied materials. Platinum or Gold (60 nm)
allowed 1 min between the application of successive
is deposited using DC magnetron sputtering. TiN
voltage pulses to allow the electrode to equilibrate,
(60 nm) is reactively sputtered from a Ti target using
which we verified by applying the voltage pulse of the
two IBS, the latter generating a nitrogen ion beam,
same voltage and making sure that the current tran-
according to [37]. IrOx was obtained via DC reactive
sients perfectly overlap. Typical electric stimulation
magnetron sputtering in an Ar/O2 plasma (100 nm)
protocols involve voltage or current pulses; therefore,
according to previous published methods [38]. ITO
the resulting current decay traces are representative
on glass was used as purchased, from Kintec, 10
for the common stimulation scenarios and were used
Ohm/square. In short, all materials consisted of a
later for model adjustment and validation.
thin film deposited by a PVD method of choice on
FA measurements were performed using a
25.4 × 25.4 mm2 cleaned glass substrates. Active elec-
PalmSens4 potentiostat and a digital oscilloscope
trochemical electrode area was defined by masking
or an Ivium PocketStat2 potentiostat. Both mod-
the substrates with adhesive polyimide stencils. The
els of used potentiostats exhibit an operational fea-
electrode area on the stencil was defined by laser cut-
ture manifested as a short delay in transient current
ting a circle with the diameter of 2 mm. Additionally,
measurement immediately after a voltage pulse. For
for validation of extrapolation experiments stencils
PalmSense4 that delay is approximately 110 µs and
with diameters of 0.5, 1, 2, 4 and 6 mm were applied
for PocketStat2 it is around 165 µs. Therefore, when
on prepared electrodes.
using PalmSens4 the current transients were meas-
ured by a digital oscilloscope, measuring a voltage
2.3. Electrochemical measurement setup drop on a 100 Ohm resistor placed in series with the
Electrodes were placed in one side of a double- potentiostat and the electrochemical cell, while with
tank electrochemical cell (MM double-tank cell, the PocketStat2 with a longer pre-trigger period the
Redox.me). The working electrode (WE) with an act- delay could be avoided.
ive area of 3.25 mm2 was fixed on one side, while
the counter electrode (CE) was placed on the other 2.6. Cyclic voltammetry (CV)
side, separated by 20 mm. The electrolyte used for CV scans were performed before the EIS measure-
the measurements (1× PBS) was prepared from ments to estimate the capacitive window for each
Roti-CELL 10x Dulbecco’s Phosphate-Buffered Saline material and after the EIS measurements to ensure
without Ca and Mg by diluting it in 1:9 ratio with DI that there were no significant changes in electrode
water according to the manufacturer specifications, performance. In addition to that, in the case of
giving a 140 mM concentration of Cl- in the final iridium oxide, repeated CV scans were performed

4
J. Neural Eng. 21 (2024) 046049 A Opančar et al

to electrochemically activate the electrode until we where the parameters Q, RS , RCT and α are obtained
could see a stable cyclic voltammogram. CV scans by fitting the measured EIS data to the Randles circuit
were performed with Palmsense4 potentiostat from (figure 2(b)).
−0.9–0.8 V vs Ag/AgCl reference wire with a scan rate If we assume that the electrode surface is homo-
of 1 V s−1 . geneous i.e. that the local surface properties are
identical for the entire electrode area, then each
2.7. Steady-state parameter extraction microscopic surface element has the same local equi-
Measurement of the electrochemical impedance spec- valent circuit describing it as in figure 1(c). Since all of
tra enables the estimation of the equivalent circuit the local equivalent circuits are connected in parallel,
parameters for different electrodes, by numerical fit- the admittance of the entire electrode corresponds to
ting of the measurements to a selected equivalent the integral of the local circuits’ admittance across the
circuit model. Over time, the EIS instrumentation entire electrode area.
has evolved to the state where portable EIS-enabled ˆ ˆ
potentiostats became a common laboratory com- Ytotal = Ys dA = Ys dA = Ys A,
modity, with the included software greatly simplify-
ing the experiments and data evaluation. Thus, our where the second equality is due to the homogeneity
focus is on using simple EIS measurements to model of the electrode.
and obtain equivalent circuit parameters for selec- Thus, we can calculate the local electrode
ted electrodes using the appropriate equivalent cir- properties for both microscopic electrode descrip-
cuit. To find an equivalent circuit that would accur- tions, Randles CPE (QSPEC , σCT , α) and Voigt RC
ately predict the electrode behaviour for a wide range (CSPEC , σCT ) by dividing the total electrode admit-
of voltages and frequencies, voltage-biased EIS meas- tance by the area of the electrode, obtaining the local
urements were performed on all the samples. electrode parameters:
In the simplest approximation, the entire elec-
trode surface is modelled as one equivalent cir- Q CDL RCT −1
cuit. However, realistic electrodes may not have QSPEC = , CSPEC = , σCT = .
A A A
an equipotential surface or may not be spatially
homogeneous—thus, it is advantageous to represent It is important to note that the series electro-
the electrode as a distribution of equivalent circuits. lyte resistance RS emerges from the geometry of
The former approach will be used in modelling the the FEM simulation and the electrolyte conductiv-
EIS data, while we will use the latter approach in our ity κ and is not implemented as an electrode prop-
FEM numerical model. erty. For actual simulations, we add an additional
Firstly, we modelled the measured data for 50Ω series resistance to compensate for the contact
each voltage bias VP to a simple RC Voigt circuit resistance of the electrode and 100Ω series resist-
(figure 2(a)). A linear RC Voigt circuit is simple ance that represents shunt resistor used in the FA
to implement in numerical simulations, and thus is measuring setup. We choose the added resistance to
widely used. However, it does not consider the reac- replicate high-frequency series resistance measured
tion kinetics, which cannot be modelled by a simple by EIS.
RC circuit, and may introduce significant errors in
a numerical simulation. However, it may be appro- 2.8. COMSOL model implementation
priate to model some of the electrode materials in a The WE, CE and REF, as well as the surrounding elec-
limited range of applied voltage. Secondly, we mod- trolyte were implemented in COMSOL Multiphysics
elled the measured data to a Randles circuit with a 6.0, using the Electric Currents interface of the
CPE (figure 2 (b)). The Randles circuit with the CPE AC/DC module. We exploited the symmetry of the
considers the reaction kinetics and diffusion beha- problem, solving it in 2D axisymmetric geometry,
viour and has provided a much better fit for all mater- which greatly reduces the complexity of the calcu-
ials. From the best numerical fit, we obtained the val- lated problem in comparison to the full 3D model.
ues of RCT (VP ), Q (VP ) and α (VP ) , as well as of the WE and CE were implemented as the 2D boundar-
electrolyte resistance Rs . From the obtained voltage ies with an appropriate boundary condition distrib-
dependent equivalent circuit parameters RCT (VP ), uted along the electrode surface, while the electrolyte
Q (VP ) and α (VP ) we calculate the distributed elec- was modelled as a 3D isotropic conductive medium
trode parameters (figures 1(b) and (c)). To obtain the characterized by the relative permittivity ϵr = 77 and
Voigt RC parameter from the single Randles equival- conductivity κ = 1.5 S m−1 or 0.15 S m−1 in the case
ent circuit, we follow the approach outlined by Brug of 1x PBS and 0.1x PBS respectively. The REF elec-
et al [30, 39] and use the formula: trode is defined as an equipotential region where the
 (1−α)/α electric VREF potential is probed, connected to a very
1 RS RCT high impedance potentiostat input. The complete
CDL = Q α ,
RS + RCT geometry with the boundary conditions is shown

5
J. Neural Eng. 21 (2024) 046049 A Opančar et al

Figure 1. A schematic workflow of methodology for modelling an arbitrary electrode and evaluating the model. Firstly, the
electrode is characterized by electrochemical impedance spectroscopy (EIS), with different constant potential (DC) biases
extending beyond the limits of the electrochemical window for a given material. In addition, current response to constant voltage
pulses for each of the materials is measured by fast amperometry (FA). EIS parameters are fitted to a chosen equivalent circuit
electrode representation, with the results implemented as a distributed electrode interface in the 3D FEM software (COMSOL
Multiphysics), in the time domain. The current response to a voltage pulse is evaluated from the FEM model, and compared to
the FA measurements. The EIS parameters of the model are iteratively optimized to obtain the best fit of a model to the FA
measurements. (a) EIS measurements for different DC biases, with Z representing the measured complex impedance (b) fitting
of the measured EIS data to the appropriate equivalent circuit, with RCT representing the charge transfer resistance, RS the
electrolyte series resistance, CPE a constant phase element with the CPE parameter Q and the CPE exponent α; (c) calculating the
distributed (intensive) electrode parameters for the Randles equivalent circuit electrode representation (CPE parameter QSPEC
and the CPE exponent α), or Voigt equivalent circuit CSPEC , and the charge transfer conductance σCT from the single equivalent
circuit (extensive) parameters, (d) implementing the distributed parameters in the 3D FEM time-dependent COMSOL model
with the control voltage VP delivered by a potentiostat between the working electrode (WE) and the reference electrode REF, by
adjusting the bias between the WE and the counter electrode CE; (e) simulating the current response (CPE FIT and CPE EIS) to
the short voltage pulse and comparing to the measurements obtained by Fast amperometry.

Figure 2. Schematic representations of (a) Voigt and (b) simplified Randles equivalent circuits, with CDL representing
electrochemical double-layer capacitance, RCT charge transfer resistance, RS electrolyte series resistance, CPE a constant phase
element with CPE parameter Q and the CPE exponent α.

6
J. Neural Eng. 21 (2024) 046049 A Opančar et al

Figure 3. 2D and the corresponding 3D rotated axisymmetric geometry with mesh depicted in the 2D image and denoted
boundary conditions and governing equations in the 3D. Current density is visualized by red arrows and lines.

in figure 3, showing the 2D and the correspond- each mesh element adjacent to the electrode surface
ing 3D rotated axisymmetric geometry. The volume whereas VM is approximated as equipotential for the
was meshed in COMSOL Multiphysics using the entire electrode. This is justified because the conduct-
Delaunay Triangular mesh with element size ranging ivity of the electrode metal is much greater than the
from 0.05 mm on the WE to 0.5 mm on the CE. conductivity of the electrolyte so the difference in
The Robin boundary condition on the WE is electric potential on the electrode side can be neg-
defined by the local normal current density on the lected compared to the potential variation along the
WE JWE [40]: electrode on the electrolyte side of the interface.
In case of the Randles circuit, σWE corresponds to
n̂·J = JWE . the surface charge density of the CPE σCPE which is
given by the expression:
The expression for JWE is calculated for every sur-
face mesh element of the WE and has two contribu- ˆt
QSPEC VWE
tions, capacitive charge density accumulating on the σWE = σCPE = dτ
Γ (1 − α) (t − τ )α
electrode surface σWE and the Faradaic current dens- 0
ity transferred across the electrode interface JF : t > 0; 0 < α < 1 , (1)
dσWE where QSPEC and α are CPE parameter and exponent,
JWE = + JF .
dt Γ is the gamma function and t is the time from the
JF and σWE are calculated according to the specif- beginning of the simulation [41]. QSPEC and α may
ics of the equivalent electrical circuit used. For the be functions of VWE . What makes this implement-
Voigt RC circuit, σWE is the surface charge density ation much more complicated and computationally
σC of the ideal capacitor given by: demanding than the Voigt circuit implementation is
the explicit dependence of the integrand on time in
σWE = σC = CSPEC VWE , the above expression. The explicit time dependence
forces us to solve the above integral anew for every
where VWE is the voltage on the electrode-electrolyte mesh element in the electrode for each time step
interface along the WE and CSPEC is the specific capa- of the solver, which makes the computation much
citance of the electrode which may be a function of slower and gives rise to additional constraints on the
VWE . time stepping while solving as well as some diffi-
We calculate VWE as the difference of the elec- culties specific to COMSOL outlined in section 2.9
tric potentials across the electrode-electrolyte inter- and the supporting information sections 1 and 2. A
face VWE ≡ VM − Vel , where Vel is the electric poten- possible simplification to significantly speed up the
tial of the electrolyte mesh elements adjacent to the computation would be considering only a single CPE
WE and VM is the potential of the metal electrode, element that represents the entire electrode instead of
which is in our case grounded via the potentiostat the distributed CPE mesh. That would require solving
circuit. Vel is spatially dependent and calculated for the integral in equation (1) only once per time step

7
J. Neural Eng. 21 (2024) 046049 A Opančar et al

instead of solving it for every electrode mesh element. of the computed potential available at each time
On the other hand, that kind of simplification would step for every mesh element of the electrode. We
not allow us to get the accurate electrode response at achieve this by mapping the entire computed his-
frequencies higher than flim or equivalently times tory VWE (r, t) into a spatial variable U (r, z, t) that we
shorter than 1/flim as demonstrated in section can then integrate over with spatial integration meth-
5.1 [35, 42]. ods. For this mapping purpose we added the stabil-
In both cases, the Faradaic current density is given ized convection–diffusion equation interface of the
by the Ohms law: COMSOL Mathematics module and constructed a
dedicated geometry for potential mapping, depicted
JF = σCT VWE, in the supporting information figure S1.
We define the 1D convection equation for the
where σCT is the charge transfer conductivity, a func-
mapping variable U:
tion of WE voltage VWE .
Since the reference electrode is connected to a very ∂U (r, z, t) ∂U (r, z, t)
high impedance input of a potentiostat we require +β =0
∂t ∂z
that there is no net current flowing through the ref-
erence electrode surface ∂REF: and we set the value of convection coefficient
β to −1 m s−1 . This means that we will have
˛
1:1 temporal to spatial mapping where U (r, z, t) =
n̂·JdA = 0.
U (r, z + βτ, t + τ ).
∂ REF We connect the variable U to the potential VWE
On the counter electrode such potential VCE is with boundary and initial conditions. We define a
applied that the voltage between WE and REF is pre- Dirichlet boundary condition:
cisely the control voltage set by the potentiostat VP .
U (r, 0, t) = VWE (r, t)
The rest of the boundaries have the Electric insulation
boundary condition given by: and the initial condition:
b
n · J = 0. U (r, z, 0) = VWE (r, 0) = 0
The equation for the electric potential in the elec- because we start the simulations with zero voltage on
trolyte bulk is derived by taking the gradient of the the electrode.
Ampere’s law and considering the identity ∇ · ∇ × The connection of U to VWE is then given by:
B ≡ 0:
  WWE (r, τ ) = U (r, 0, τ ) = U (r, z = β (t − τ ) , t) .
∂D
∇ Jc + = 0.
∂t Finally, by substitution we can rephrase the integ-
ral in equation (1) in terms of spatial integration of
After relating the conduction current density in
variable U:
the electrolyte Jc to electrical field with Ohm’s law
with bulk conductivity κ: ˆβ t
QSPEC U (r, z, t)
σCPE (r, t) = dz t > 0; 0 < α < 1.
Γ (1 − α) β(t − z/β)α
Jc = κE 0

and taking the constitutive relation for D: We evaluate this integral using the linproj operator
in COMSOL 6.0. Additional details of the implement-
D = ϵ0 ϵr E,
ation can be found in the supporting information
finally, the electric field is connected to the electric sections 1 and 2. The COMSOL model file, pdf report
potential by a negative gradient: and supporting data can be obtained from the repos-
itory[43] referenced in the data availability statement.
E = −∇V.
2.10. Transient parameter optimization
2.9. Time-domain implementation of the CPE In addition to EIS, we have also measured the cur-
The additional difficulty with the CPE implementa- rent response to constant voltage pulses applied to
tion arises from the fact that coupling of the com- the electrodes by FA and by using those measure-
puted potentials at former time steps back into cal- ments we evaluate our final electrode model. We can-
culation at current time step via integration is not not expect that the equivalent circuit with paramet-
possible in COMSOL 6.0. to the best of our know- ers obtained from EIS will work perfectly for the fast
ledge. That makes the implementation of integral voltage pulse transients in the time domain, as the
in equation (1) impossible with available methods biased EIS measurement protocol establishes the DC
for time integration within COMSOL. To circumvent equilibrium state before starting AC measurement,
that issue, we needed to make the entire history and thus it may miss some of the dynamic processes

8
J. Neural Eng. 21 (2024) 046049 A Opančar et al

that play a key role in the case of fast pulses, such as the simple and numerically robust. Voltage dependent
dynamics of the oxygen reduction and other Faradaic electrode parameters QSPEC (VWE ) , CSPEC (VWE ) and
reactions at the electrodes. The reason for that is that σCT (VWE ) are implemented in COMSOL by using the
EIS measurements can, for example, quickly deplete built-in interpolation function feature which creates
oxygen locally available at the electrode even before an analytic function from the parameter values meas-
the impedance measurement begins, and as a result ured or fitted with some finite voltage resolution.
the process becomes diffusion limited, which is not These interpolation functions are used in the physics
the case to the same extent for the voltage pulses at governing equations as the parameter values. In this
the millisecond time scale [18]. We can, however, find way the parameter values are dynamically updated at
the equivalent circuit which will accurately represent every solver time step to have the value prescribed by
the electrode in the stationary EIS limit and hope- the functional dependence given the value of VWE at
fully, by small adjustments, we can have a good fit on that exact time step.
the current transients measured by the FA. Thus, to
improve the model’s results, we did a series of numer-
ical optimizations of the equivalent circuit paramet- 3. Results and discussion
ers starting from the EIS data, by using BOBYQA
algorithm in the Parameter Estimation study of the We have built our 3D FEM model to faithfully rep-
Optimization module of COMSOL 6.0. We varied the resent the geometry of the measured electrochem-
values of QSPEC and CSPEC to obtain the best least- ical cell. The fundamental microscopic (material-
squares fit to the measured data. dependent) electrode parameters were parametrized
This adjustment is not trivial because the circuit by two equivalent circuits, Voigt (RC) and Randles
parameters QSPEC (VWE ) and CSPEC (VWE ) are func- (CPE), applied to the electrode mesh elements. While
tions of electrode voltage and for each pulse the the choice of the Voigt RC circuit may be naïve, it is
voltage on the electrode changes from VWE (t = 0) very commonly used due to the ease of implement-
= 0 to VWE (t = tP ) ≈ VP where tP and VP refer to the ation, and we will show in which cases it may still
pulse duration and pulse voltage respectively. That be used. On the other hand, the Randles equivalent
means that for each voltage pulse we need to fit a circuit can model the measured electrodes faithfully,
function rather than a simple number for QSPEC and while unfortunately introducing increased complex-
CSPEC which drastically complicates the fitting pro- ity to the model’s implementation. There are a few
cedure and does not guarantee a unique result for key insights from these results that will guide us in
each pulse due to possible local minimums that the the choice of the most appropriate equivalent circuit.
fitting procedure may encounter. To solve this issue, First, in all of the EIS measurements, we see the pres-
we use the following fitting procedure. ence of the CPE characterized by ‘depressed semi-
First, we fit to the voltage pulses of the smallest circles’ in the Nyquist plot or flat plateaus visible in
amplitude VP = VP±1 = ±100 mV vs Ag/AgCl elec- the phase component of the Bode plot shown in the
trode for anodic and cathodic pulse. We can fit the supporting information figures S2–S7. That implies
values for QSPEC and CSPEC for the limiting VWE = that we must use an equivalent circuit with a CPE ele-
VP±1 only and assume that for the voltages VP−1 < ment, such as a Randles circuit, if we want to obtain
VWE < VP+1 the parameters QSPEC and CSPEC change an accurate fit. Secondly, the values of the fitted para-
linearly between the end values. That assumption meters for the same material vary greatly for different
is well justified if the electrode parameters do not DC biases. That means that we must allow for voltage-
change drastically within the voltage resolution of our dependent parameters in our equivalent circuit.
measurements ∆V = |VP+1 − VP−1 |. As mentioned in 2.1, the observed CPE behaviour
Next, we fit on the pulses of the second smal- can be the result of a 2D or 3D variation of electrode
lest amplitudeVP = VP±2 = ±200 mV vs Ag/AgCl. properties or, above the limiting frequency flim , the
We again fit the values for QSPEC and CSPEC for the result of the inhomogeneous electric field distribu-
limiting values VWE = VP±2 only and assume that tion along the electrode surface. If the reason for the
the parameters change piecewise linearly between observed CPE behaviour is one of the former two,
VP−2 < VWE < VP−1 and VP+1 < VWE < VP+2 where there is no other choice than to use the Randles circuit
the values for the parameter between VP−1 and VP+1 (figure 2(b)) both for fitting the EIS data (figure 1(b))
are fixed from the previous two pulses and are not and as a distributed microscopic description of the
subject to further change. electrode surface that we implement in the FEM
We thus repeat the same procedure for all the simulation (figure 1(c)). If the latter is the reason
voltage pulses in the ascending order of voltage for the observed CPE behaviour, the Randles circuit
amplitude. In this way we only fit a single value must still be used to fit the EIS measurements, but
for QSPEC and CSPEC for each pulse and reduce there is the possibility that the microscopic electrode
the fitting procedure to a series of single variable parameters can be well described by a simpler Voigt
fits, making the entire procedure computationally RC circuit (figure 2(a)) and that we will obtain the

9
J. Neural Eng. 21 (2024) 046049 A Opančar et al

Figure 4. A Nyquist plot of the measured, fitted and FEM simulated EIS spectrum for a gold electrode, demonstrating the
advantage of a distributed equivalent circuit in a FEM model in comparison to using a single domain EIS equivalent circuit in the
high-frequency region.

observed CPE behavior ‘for free’ as a result of the real- measured spectrum was found for Au, as well as for
istic 3D FEM distributed electrode model instead of other materials, with the high frequency region for the
the single equivalent circuit. Au electrode shown in the inset. We can see that all
We will not try to a priori determine if the micro- the data are in good agreement for lower frequencies,
scopic RC description is completely justified or not; tracing almost a perfect straight line in the Nyquist
rather we will present both implementations and plot characterized by the CPE coefficient α = 0.921.
report on the relative error of using each micro- Above the limiting frequency (426 Hz for the gold
scopic description. The reader then has to determ- electrode) the measured and simulated distributed
ine whether the relative error associated with each data deviates from the trend, with significant differ-
description is acceptable for their specific application. ences above three times the limiting frequency. The
Once our model is verified by comparing it to reason is that a 3D model can capture the variations in
the measurements, the electrode parameters can be potential close to the electrode edge, whereas a single
applied to arbitrary electrode geometries, such as in equivalent circuit cannot. In our FEM model every
multi-electrode arrays or other electrodes of choice. mesh element behaves as a separate equivalent circuit;
We will demonstrate the power of extrapolation of the thus, our model is shown to reproduce the physics
model by applying it to different geometries, stimula- more realistically than a simple single domain equi-
tion protocols and electrolyte environments. valent circuit obtained from the EIS [35, 42].

3.1. Model verification by EIS 3.2. Model verification by FA

To verify our implementation of the 3D electrode Electrical stimulation protocols commonly employ
model in the time domain, we firstly replicate the pulsed stimulation, which displays complicated tem-
results obtained from EIS by doing the time domain poral dynamics changing from kinetic to diffusion
simulations for the small AC perturbation and using controlled regimes as the pulse progresses—in con-
the inverse Fourier transform to get the frequency trast to the EIS measurements, which show responses
dependent impedance data. We calculate the local to single sinusoidal frequencies after the DC equilib-
equivalent circuit parameters (section 2.7) from fit- rium is established. Thus, the model needs to be veri-
ting the EIS data to single Randles circuit (figure 2(b)) fied for the intended usage scenario. Voltage pulses
and implement them in our model in a distributed within a range of amplitudes, typically within the
form (figure 1(c)). We compare the numerically mod- electrochemical window of water-based electrolytes,
elled result with the impedance spectra fitted to a were applied to all electrodes, and current transients
single Randles circuit and the measured data. Here were recorded.
we show only the representative results for a gold Firstly, we used the voltage dependent paramet-
electrode (figure 4). A good match over the whole ers for the distributed Randles and Voigt equivalent

10
J. Neural Eng. 21 (2024) 046049 A Opančar et al

Figure 5. Simulated and measured current traces as a response to a 5 ms square −900 mV voltage pulse for the TiN electrode. The
current trace modelled from the parameters obtained by EIS is underestimating the transferred charge. Deviations from the
measured data are shown in the upper panel.

circuits obtained from EIS and plugged them dir- and Voigt circuits that provide the best fit on the FA
ectly into our model. As can be seen in figure 5, pulses on each of the measured materials and report
the resulting traces, even for the CPE fit, differed the errors for the two implementations.
significantly from the measured data, most likely EIS measurements can provide a good paramet-
due to the dynamics of the electrode which was rization for different electrodes operated in pulsed
not captured fully by the EIS measurements. To mode near their open circuit potentials (OCP). For
improve the model’s results, we did a series of numer- pulse voltages further from their OCPs the EIS para-
ical optimizations of the equivalent circuit paramet- meters can serve as a good starting point for further
ers starting from the EIS data. The comparison of optimization.
voltage-dependent Voigt circuit parameters (σCT and
CSPEC ) and Randles circuit parameters (QSPEC and α) 3.3. Voigt versus Randles implementation
obtained from EIS, and optimized to the measured As mentioned, the Voigt RC equivalent circuit is often
data, starting from the EIS parameters, are shown in implemented in time-domain simulations, due to the
figure 6. rather simple numerical implementation.
We observe large differences between steady state To quantify the agreement between measured and
parameter obtained by EIS and transient paramet- simulated data we introduce the error metric as a
ers obtained by fitting to FA measurements in all measure of relative difference between the simulated
the measured materials at some potential regions. We and measured delivered electric charge:
hypothesize that in all cases the difference originates ∫ |IS − IM | dt
from Faradaic reactions which will not be apparent in Error =
∫ |IM | dt
EIS to the extent which they are visible in FA because
in EIS method the system must assume a steady state where IS and IM are simulated and measured current.
at given DC bias before applying AC perturbation. We calculated the error metric for each voltage pulse
This DC equilibrium will deplete reactants close to for every material (figure S9). The error metric aver-
the electrode thus making a reaction seem much less aged across voltages for different materials is shown in
prominent than in a short pulse transient where the figure 7 with error bars indicating standard deviation
reactants are present in the equilibrium amount near of the error metric for the different voltages.
the electrode surface. In particular, gold electrode in With the Randles circuit (CPE) it is possible to
the cathodic region showed almost five times greater obtain good agreement between measured and simu-
CPE parameter in transient than predicted by EIS lated data with an average error below 10% for all the
and we can attribute that difference specifically to the measured materials. With the Voigt circuit (RC), the
to the oxygen reduction reactions according to [18]. average error is seen clearly increasing with decreas-
Finally, we present parametrizations for the Randles ing CPE exponent α. Great care should be taken when

11
J. Neural Eng. 21 (2024) 046049 A Opančar et al

Figure 6. (a) Results of best fits of the voltage dependent (a) RC parameters—specific double layer capacitance CSPEC and charge
transfer conductivity σCT , and (b) CPE parameters QSPEC and αCPE , for titanium nitride, titanium, indium tin oxide, gold,
platinum and iridium oxide. Parameters obtained from EIS are labelled EIS, fitted parameters that offer improved agreement with
the measured data are labelled FIT. Voltages were measured vs. Ag/AgCl electrode. Error bars indicate fit uncertainty.

12
J. Neural Eng. 21 (2024) 046049 A Opančar et al

Figure 7. Average error (∆) of modelled in comparison to the measured current traces for the RC and the CPE mode, for
different materials. Error bars indicate standard deviation for all the voltage amplitudes for each material.

applying RC parametrization to materials with α ⩽ that the model works well for the measured variation
0.9 as the possible error is above 25%. in electrode sizes that spans 12 times ratio in diameter
or 144 times ratio in surface area between the largest
3.4. General model applicability and the smallest measured electrode size.
To determine the applicability of our model in more We do observe a significant difference between
general cases as well as the limitations of the model the agreement for different electrolyte concentrations
we performed additional experiments varying elec- with 1xPBS, in which the parameter extraction was
trode size, electrolyte concentration and stimulus performed, outperforming 0.1x PBS in average error
shape and thus tested its extrapolation validity for cer- by about a factor of two. Even though that difference
tain cases. We chose to focus these experiments on a is significant, it can be put into perspective of the rel-
TiN electrode as a representative material due to the ative errors of different electrode model implement-
widespread use for electrical stimulation as well as its ations (Voigt RC vs. CPE) where the RC model still
median position for CPE exponent (figure 7). performs significantly worse (figure 7) than the CPE
model for a 0.1x PBS electrolyte concentration.
3.4.1. Variations in electrode size and electrolyte It is apparent in figures 8(b) and (c) that for
concentration 0.1x PBS the model systematically overestimates the
We performed FA measurements for the combina- current indicating the overestimation of the double
tions of 5 different electrode sizes (with diameters 0.5, layer capacitance. The same is also be predicted
1, 2, 4 and 6 mm) and two electrolyte concentrations from theory of the Gouy–Chapman–Stern model of
(1x PBS and 0.1x PBS). The parameters are extrac- the double layer capacitance where the capacitance
ted and fitted as explained in Chapter 4 only for the increases with electrolyte concentration but levels off
2 mm electrode in 1x PBS. These parameters are then at high electrolyte concentrations [44]. Since we are
applied to all the other size/electrolyte combinations, doing 10x dilution of the electrolyte we might be see-
simulated in COMSOL, and compared with the FA ing this reduction in double layer capacitance which
measurements. Comparison between simulated and the model is incapable of predicting since we gave it
measured currents as a response to a short voltage the parameters extracted from the measurement in 1x
pulse is shown in figure 8. PBS.
We do not observe a particular trend in the agree- Although it might be possible to explicitly imple-
ment between measured and simulated currents for ment the compensation of the double layer capa-
the different electrode diameters for neither of the citance change with electrolyte concentration in the
electrolyte concentrations as seen in figure 8(d)). We model from the theory, we do not think that would
also do not see the combination on which the para- be very worthwhile for two reasons. First, the relat-
meter extraction was performed as particularly stand- ive error for 0.1x PBS is still not that large (com-
ing out from the rest of the data so we can conclude pared to Voigt RC implementation for example) and

13
J. Neural Eng. 21 (2024) 046049 A Opančar et al

Figure 8. Comparison between simulated and measured currents with a TiN electrode as a response to a short voltage pulse for
different electrode diameter and electrolyte concentration combinations. (a) Measured and simulated currents for different
electrode diameters in 1x PBS and voltages 0.8 V (positive) and −0.9 V (negative) vs Ag/AgCl reference. Parameters are extracted
from the 2 mm electrode (shown in red) and applied to all the other combinations. (b) Same as in (a) but for a 0.1× PBS. (c)
Positive current traces from figures (a) and (b) combined in a Log–Log plot to provide improved visibility with largely differing
current and time scales. (d) Average error as defined in 5.3 for different diameter-electrolyte combinations. Error for each
combination is averaged for 17 voltage pulse amplitudes from −900 to 800 mV in steps of 100 mV vs Ag/AgCl. The error bar
represents standard deviation for the 17 voltage amplitudes for each combination.

the model can still be considered usable for elec- of current amplitudes and adjusted the pulse duration
trolyte concentrations variations not larger than 10x such that the maximum voltage on the electrode does
which corresponds to a significant conductivity range not exceed −900 mV vs. Ag/AgCl reference which was
in which many typical electrical stimulation experi- the maximum value for all voltage-controlled exper-
ments would fall into. Second, it would be much sim- iments. The same model parameters as in the previ-
pler to just use the electrolyte preparation for elec- ous subsection are used without any further optimiz-
trode characterization that is as similar as possible to ation. Simulated and measured voltages are shown in
the one used in the in vitro or in vivo experiments. figure 9.
Overall, we observe that the measured and sim-
3.4.2. Variation in stimulus shape ulated voltages are in agreement in both the cath-
Up to now, we have only used voltage-controlled odic and the anodic part of the pulse and the agree-
pulses because it is relatively simple to determine the ment holds for a range of current amplitudes and
safe voltage range from the CV that will not cause pulse durations. The error metric as defined in 5.3
excessive Faradaicity that can for example cause sig- with current replaced by voltage gives us an average
nificant formation of bubbles on the electrode which error of (7 ± 2) % which is comparable to the error
would make modelling very difficult or outright dam- for voltage pules (figures 7 and 8). We can therefore
age the electrode. In practice though, researchers are conclude that the model holds similarly well for cur-
often interested in the deposited charge which is rent controlled stimulation.
much simpler to control with the current controlled
pulses. 3.5. Interpretation of model parameters and model
Thus, we tested the ability of our model to rep- applicability
licate the electrode behaviour for current controlled For bioelectronics applications, be it for electrical
pulses on biphasic charge balanced cathodic leading stimulation or electrical measurements, electrode
pules as a representative waveform. We chose a range impedance is of essential value. As a rule of thumb,

14
J. Neural Eng. 21 (2024) 046049 A Opančar et al

Figure 9. Simulated and measured voltages for a 2 mm diameter TiN electrode in current controlled pulses for different current
amplitudes and pulse durations.

low impedance is desirable in most cases, as is high 4. Conclusions

interfacial capacitance and low Faradaic current com-
ponent. IrOx and TiN emerge as good electrode We outlined a method that can be used to estim-
candidates in this respect. However, subtle effects ate and simulate electrode response at the elec-
of electrode geometry and cell-electrode cleft dis- trode/electrolyte interface for pulsed operation from
tance may have a profound effect on the stimula- electrode parameters that can be easily obtained by
tion thresholds, as well as on measured signal amp- the common experimental techniques of EIS and FA
litudes. This work presents a tool to be used in com- .This method gives a sub-10% error between meas-
prehensive numerical models, which would include ured and simulated currents for a range of materials
both the presented realistic electrode model, and a (Ti,TiN,ITO,Au,Pt,IrOx) and pulsed voltages when
model of a biological system, be it a single cell, a the CPE model is used. When using the ideal capa-
neural bundle or a portion of the central nervous citor instead, the agreement is good only for materi-
system. In addition, we give an important caveat als that have a CPE exponent α greater than 0.9. On
about using the Voigt RC equivalent circuit to rep- the other hand, using the CPE has many disadvant-
resent most of the presented materials—even though ages for example, less intuitive physical interpreta-
the time-domain implementation of an equivalent tion of the electrode parameters, complicated imple-
RC circuit is simple, the absolute error in modelled mentation in time-dependent simulations, and it is
delivered charge can easily reach up to 50% for cer- computationally much more demanding than simple
tain materials. If the modelling outcomes differ signi- Voigt RC electrode representation. It is thus helpful
ficantly with the variation of delivered charge of this to understand and estimate a possible error that one
scale, a more complex model including the distrib- would make in choosing an RC representation for a
uted CPEs might be necessary. The model has shown particular electrode material and weigh the pros and
to be robust for current-driven stimulation, as well cons of choosing to use the RC or CPE model.
as when used in up to 10x diluted electrolytes. The
RC and CPE parameters we give for different mater- Data availability statement
ials should be used only as first approximations for
the electrode used, and always the fresh model should The data that support the findings of this study are
be built for every new application and every new openly available at the following URL/DOI: https://
material. urn.nsk.hr/urn:nbn:hr:217:393738.

15
J. Neural Eng. 21 (2024) 046049 A Opančar et al

Author contributions References

Aleksandar Opančar: investigation (lead); software [1] Maynard E M, Nordhausen C T and Normann R A 1997 The
utah intracortical electrode array: a recording structure for
(lead); methodology (equal); formal analysis (lead);
potential brain-computer interfaces Electroencephalogr. Clin.
visualization (lead); writing—original draft (sup- Neurophysiol. 102 228–39
porting); writing—review and editing (equal). Eric [2] Zeck G and Fromherz P 2001 Noninvasive neuroelectronic
Daniel Głowacki: conceptualization (supporting); interfacing with synaptically connected snail neurons
immobilized on a semiconductor chip Proc. Natl Acad. Sci.
writing—review and editing (equal); funding acquis-
98 10457–62
ition (supporting); resources (supporting);Vedran [3] Fromherz P 2002 Electrical interfacing of nerve cells and
Ðerek: conceptualization (lead); supervision (lead); semiconductor chips ChemPhysChem 3 276–84
writing—original draft (lead); writing—review and [4] Dinyari R, Loudin J D, Huie P, Palanker D and Peumans P
2009 A curvable silicon retinal implant Technical Digest Int.
editing (equal); methodology (equal); data curation;
Electron Devices Meeting, IEDM (Baltimore, MD, USA, 7–9
funding acquisition (lead); project administration. December 2009) (https://doi.org/10.1109/IEDM.2009.
5424291)
Conflict of interest [5] Schmidt T et al 2022 Light stimulation of neurons on
organic photocapacitors induces action potentials with
millisecond precision Adv. Mater. Technol. 7 2101159
The authors have no conflicts to disclose. [6] Silverå Ejneby M et al 2021 Chronic electrical stimulation of
peripheral nerves via deep-red light transduced by an
implanted organic photocapacitor Nat. Biomed. Eng.
Funding statement 6 741–53
[7] Silverå Ejneby M, Migliaccio L, Gicevičius M, Ðerek V,
This work has been supported by the Croatian Science Jakešová M, Elinder F and Głowacki E D 2020 Extracellular
Foundation under the Project UIP-2019-04-1753. photovoltage clamp using conducting polymer-modified
organic photocapacitors Adv. Mater. Technol.
V Ð and A O acknowledge the support of project
5 1900860
CeNIKS co-financed by the Croatian Government [8] E K 2014 Implantable Bioelectronics ed E Katz (Wiley)
and the European Union through the European [9] Deshmukh A et al 2020 Fully implantable neural recording
Regional Development Fund—Competitiveness and stimulation interfaces: peripheral nerve interface
applications J. Neurosci. Methods 333 108562
and Cohesion Operational Programme (Grant
[10] Won S M, Cai L, Gutruf P and Rogers J A 2023 Wireless and
No. KK.01.1.1.02.0013), and the QuantiXLie battery-free technologies for neuroengineering Nat. Biomed.
Center of Excellence, a project co-financed by the Eng. 7 405–23
Croatian Government and European Union through [11] Cai L and Gutruf P 2021 Soft, wireless and subdermally
implantable recording and neuromodulation tools J. Neural
the European Regional Development Fund—
Eng. 18 41001
the Competitiveness and Cohesion Operational [12] Lanmüller H, Sauermann S, Unger E, Schnetz G, Mayr W,
Programme (Grant KK.01.1.1.01.0004). This work Bijak M, Rafolt D and Girsch W 1999 Battery-powered
was supported by the European Research Council implantable nerve stimulator for chronic activation of two
skeletal muscles using multichannel techniques Artif. Organs
(ERC) under the European Union’s Horizon 2020
23 399–402
research and innovation program (E.D.G. Grant [13] Cogan S F, Ludwig K A, Welle C G and Takmakov P 2016
Agreement No. 949191), by the Grant Agency Tissue damage thresholds during therapeutic electrical
of the Czech Republic under Contract No. 23- stimulation J. Neural Eng. 13 021001
[14] Shannon R V 1992 A model of safe levels for electrical
07432S, and by funding from the National Center
stimulation IEEE Trans. Biomed. Eng. 39 424–6
for Neurological Research, supported by the Czech [15] Reilly J P 1998 Applied Bioelectricity (Springer)
Ministry of Education, Youth, and Sports (MEYS CR, [16] Merrill D R, Bikson M and Jefferys J G R 2005 Electrical
LX22NPO5107). Sample fabrication was enabled by stimulation of excitable tissue: design of efficacious and safe
protocols J. Neurosci. Methods 141 171–98
CzechNanoLab Research Infrastructure supported by
[17] Merrill D R 2010 The electrochemistry of charge injection at
MEYS CR (LM2023051). the electrode/tissue interface Implantable Neural Prostheses 2
ed D Zhou and E Greenbaum (Springer) pp 85–138
Ethics statement [18] Ehlich J, Migliaccio L, Sahalianov I, Nikić M, Brodský J,
Gablech I, Vu X T, Ingebrandt S and Głowacki E D 2022
Direct measurement of oxygen reduction reactions at
Ethics approval not required. neurostimulation electrodes J. Neural Eng. 19 036045
[19] Günter C, Delbeke J and Ortiz-Catalan M 2019 Safety of
ORCID iDs long-term electrical peripheral nerve stimulation: review of
the state of the art J. Neuroeng. Rehabil. 16 13
[20] Harris A R 2020 Current perspectives on the safe electrical
Aleksandar Opančar  https://orcid.org/0000-0003- stimulation of peripheral nerves with platinum electrodes
3471-1110 Bioelectron. Med. 3 37–49
Eric Daniel Głowacki  https://orcid.org/0000- [21] Díaz L et al 2020 Ethical considerations in animal research:
the principle of 3r’s Rev. Invest. Clin. 73 199–209
0002-0280-8017 [22] Singh M, Fuenmayor E, Hinchy E P, Qiao Y, Murray N and
Vedran Ðerek  https://orcid.org/0000-0001-9507- Devine D 2021 Digital twin: origin to future Appl. Syst.
6865 Innov. 4 36

16
J. Neural Eng. 21 (2024) 046049 A Opančar et al

[23] Wright L and Davidson S 2020 How to tell the difference [34] Jorcin J B, Orazem M E, Pébère N and Tribollet B 2006 CPE
between a model and a digital twin Adv. Model. Simul. Eng. analysis by local electrochemical impedance spectroscopy
Sci. 7 Electrochim. Acta 51 1473–9
[24] Zimmermann J, Budde K, Arbeiter N, Molina F, Storch A, [35] Huang V M, Vivier V, Orazem M E, Pe´bère N and Tribollet B
Uhrmacher A M and van Rienen U 2021 Using a digital twin 2007 The apparent constant-phase-element behavior of an
of an electrical stimulation device to monitor and control the ideally polarized blocking electrode J. Electrochem. Soc.
electrical stimulation of cells in vitro Front. Bioeng. 154 C81
Biotechnol. 9 765516 [36] Huang V, Vivier V, Orazem M E, Pebere N and Tribollet B
[25] Ghazavi A, Westwick D, Xu F, Wijdenes P, Syed N and 2007 The apparent cpe behavior of a disk electrode with
Dalton C 2015 Effect of planar microelectrode geometry on faradaic reactions: a global and local impedance analysis ECS
neuron stimulation: finite element modeling and Trans. 3 567–85
experimental validation of the efficient electrode shape J. [37] Gablech I, Migliaccio L, Brodský J, Havlíček M, Podešva P,
Neurosci. Methods 248 51–58 Hrdý R, Ehlich J, Gryszel M and Głowacki E D 2023
[26] Buitenweg J R, Rutten W L C and Marani E 2000 Finite High-conductivity stoichiometric titanium nitride for
element modeling of the neuro-electrode interface IEEE Eng. bioelectronics Adv. Electron. Mater. 9 2200980
Med. Biol. Mag 19 46–52 [38] van Ooyen A, Topalov G, Ganske G, Mokwa W and
[27] Zimmermann J, Sahm F, Arbeiter N, Bathel H, Song Z, Schnakenberg U 2009 Iridium oxide deposited by pulsed
Bader R, Jonitz-Heincke A and van Rienen U 2023 dc-sputtering for stimulation electrodes J. Micromech.
Experimental and numerical methods to ensure Microeng. 19 074009
comprehensible and replicable alternating current electrical [39] de Pauli M, Gomes A M C, Cavalcante R L, Serpa R B,
stimulation experiments Bioelectrochemistry 151 108395 Reis C P S, Reis F T and Sartorelli M L 2019 Capacitance
[28] Lasia A 2014 Electrochemical Impedance Spectroscopy and Its spectra extracted from EIS by a model-free generalized
Applications (Springer) phase element analysis Electrochim. Acta
[29] Lukács Z and Kristóf T 2020 A generalized model of the 320 134366
equivalent circuits in the electrochemical impedance [40] Cantrell D R, Inayat S, Taflove A, Ruoff R S and Troy J B 2008
spectroscopy Electrochim. Acta 363 137199 Incorporation of the electrode–electrolyte interface into
[30] Brug G J, van den Eeden A L G G, Sluyters-Rehbach M and finite-element models of metal microelectrodes J. Neural
Sluyters J H 1984 The analysis of electrode impedances Eng. 5 54–67
complicated by the presence of a constant phase element J. [41] Westerlund S 1991 Dead matter has memory Phys. Scr.
Electroanal. Chem. 176 275–95 43 174–9
[31] Rossmacdonald J 1984 Note on the parameterization of the [42] Newman J 1970 Frequency dispersion in capacity
constant-phase admittance element Solid State Ion. 13 147–9 measurements at a disk electrode J. Electrochem. Soc.
[32] Zoltowski P 1998 On the electrical capacitance of interfaces 117 198
exhibiting constant phase element behaviour J. Electroanal. [43] Opančar A, Ðerek V and Glowacki E 2024 Choosing the
Chem. 443 149–54 Right Electrode Representation for Modeling Real
[33] Córdoba-Torres P, Mesquita T J and Nogueira R P 2013 Bioelectronic Interfaces: A Comprehensive Guide [Data set]
Influence of geometry-induced current and potential Prirodoslovno-matematiˇcki fakultet urn:nbn:hr:217:393738
distributions on the characterization of constant-phase [44] Bard A J, Faulkner L R and White H S 2022 Electrochemical
element behavior Electrochim. Acta 87 676–85 Methods: Fundamentals and Applications (Wiley)

17