From Finite Element Simulations to Equivalent Circuit Models of Extracellular Neuronal Recording Systems Based on Planar and Mushroom Electrodes

Objective: define a new methodology to build multi-compartment lumped-elements equivalent circuit models for neuron/electrode systems. Methods: the equivalent circuit topology is derived by careful scrutiny of accurate and validated multiphysics finite-elements method (FEM) simulations that couple ion transport in the intra- and extracellular fluids, activation of the neuron membrane ion channels, and signal acquisition by the electronic readout. Results: robust and accurate circuit models are systematically derived, suited to represent the dynamics of the sensed extracellular signals over a wide range of geometrical/physical parameters (neuron and electrode sizes, electrolytic cleft thicknesses, readout input impedance, non-uniform ion channel distributions). FEM simulations point out phenomena that escape an accurate description by equivalent circuits; notably: steric effects in the thin electrolytic cleft and the impact of extracellular ion transport on the reversal potentials of the Hodgkin-Huxley neuron model. Conclusion: our multi-compartment equivalent circuits match accurately the FEM simulations. They unveil the existence of an optimum number of compartments for accurate circuit simulation. FEM simulations suggest that while steric effects are in most instances negligible, the extracellular ion transport affects the reversal potentials and consequently the recorded signal if the electrolytic cleft becomes thinner than approximately 100 nm. Significance: the proposed methodology and circuit models improve upon the existing area and point contact models. The coupling between the extracellular concentrations and reversal potential highlighted by FEM simulations emerges as a challenge for future developments in lumped-element modeling of the neuron/sensor interface.


I. INTRODUCTION
Since many years micro-electrode arrays (MEAs) [1] are the workhorse of in-vitro extracellular recording of neuronal electrical activity [2], [3].Despite their limitations compared to alternative intracellular recording methods [4], they are preferred for the longterm stability [5], CMOS-compatible fabrication process [6], and high parallelism [7], enabling measurements from large neuron populations.MEA design and the interpretation of MEA data require a detailed understanding of the recorded signal generation mechanisms.This task benefits from physics-based distributed models describing the system down to sub-micron scale in combination with reliable equivalent circuit models for the neuron/electrode interaction.
Distributed modeling of extracellular recordings typically relies on the quasistatic approximation of Maxwell equations and on the ohmic assumption for the extracellular fluid [8].Herein, the resulting system is commonly solved in a hybrid fashion: first, the currents of the neuronal membrane are computed according to the cable theory [9] assuming a grounded extracellular domain; second, they are used as sources to compute the extracellular potentials with either analytical approximations [10], [11] or finite-element methods (FEM) [12], [13].Alternatively, FEM approaches solve self-consistently the dynamics of both the intracellular and extracellular fluids with an explicit 3D geometry of the cell membrane [14]- [16].In all cases, it is possible to include a description of the sensing electrodes [13], [14].The ohmic assumption can be relaxed by resorting to an electrodiffusive description of the ionic transport phenomena occurring in the cellular fluids.This can be accomplished with hybrid approaches [17] as well as self-consistent formalisms with explicit descriptions of the neural geometries [18], [19].Electrodiffusive models allow a close inspection of the contribution of each ionic species to the recorded signal and, furthermore, they provide more accurate solutions for spatial resolutions lower than the micrometer [20], a scale which is typical of the neuron/electrode cleft [21].Notwithstanding, relatively little effort has been made in the past to model the neuron/electrode interaction at the electrodiffusive level [22].We have recently proposed a novel FEM framework [23], which sets the basis for the present work.
Equivalent circuit models are a widely used alternative to distributed models [24], [25].Indeed, the cable theory formalism mentioned above [9] is amenable for integrations in circuit simulators [26]- [28].Building upon this fact, lumped-elements equivalent circuits allowed many groups [29]- [32] to gain insight into how the biological/physical/geometrical parameters of the recording system affect the experimentally recorded signals.Most circuits rely on the point contact model [21]: the electrode and the junctional portion of the neuron membrane (i.e., that interfaced to the electrode) are lumped together into a netlist node, under the assumption that the interface is far away from the ground.The electrode is commonly described by passive elements whereas the membrane is represented with different levels of accuracy and complexity (in decreasing order): i) Hodgkin-Huxley (HH) [33] and HH-like models [29], [30], [32], [34], [35]; ii) RC passive driven by voltage/current sources with HH-like waveforms [4], [5], [36]; iii) RC passive with AC input drive [5], [31], [37].Still, the non-junctional membrane (i.e., that not facing the electrode) can be described by one or more compartments [30], [38].Ref. [35] demonstrated that the aggregation of HH blocks in a multi-compartmental description of the neuron, in combination with a mono-compartmental circuit description of the junctional neuron/electrode interface, accurately reproduces the waveforms of recorded signals.
However, the point contact model cannot account for distributed effects at the electrode/neuron junction which can modify the local extracellular potential.To this aim, area contact models [21] better replicate these effects by considering more than one compartment for the junction.These models typically describe membranes as distributed RC nets with intracellular potentials/currents left as external sources, e.g., for frequency analysis [39], or describe only the recording device in a distributed fashion while including a single HH compartment for the membrane [40].To the best of our knowledge, truly distributed area contact models with HH-based junction implementation for transient studies are still missing.Ideally, such an equivalent circuit should account for the distributed nature of the neuron/electrode junction and of the HH-type membrane, and should adopt a suitable number of compartments to preserve the accuracy of the recorded signal.To this end, physics-based FEM modeling can be used as a testbed to evaluate the accuracy of the circuit across various neuron/electrode/junction geometries and physical parameters, thus reducing the need for experimental data.Such methodology may additionally shed light on i) the minimum number of junctional and non-junctional compartments required; ii) the optimal approach to translate the physical domain into compartments.This work aims to overcome the limitations of most existing area contact models and provide a tool to bridge the gap between these points.
First, we present an extended version of the FEM model proposed in [23] to investigate continuous and non-uniform ion channel distributions, neuron membranes with domed and ellipsoidal shapes, and realistic planar and mushroom geometries of the sensing electrode.Secondly, we propose a method to derive accurate multicompartmental equivalent circuits for the geometries analyzed with FEM simulations.The circuit and FEM models are compared for various design parameters to assess the minimum number of compartments necessary for an accurate equivalent circuit description of the sensing system.Finally, we raise awareness toward situations with thin junctional clefts, where neglecting the electrodiffusive nature of the extracellular fluid fundamentally degrades the accuracy of the circuit description, due to an alteration of the reversal potentials at the neuron membrane.
The manuscript proceeds as follows.The FEM simulation framework and the derivation of the equivalent circuit are described in section II.A comparison between the two models and a systematic analysis of the influence of the system geometry on the results are reported in section III.Conclusions are finally drawn in section IV.

A. FEM modeling framework
The FEM model is an extension of the one in [23].It is based on the tool COMSOL Multiphysics ® [44] and aims to describe the intra-and extracellular fluids, the neuron membrane, and the metallic electrode coupled to a readout circuit.This is showcased in Figure 1 which also defines the model equations.All the physical variables in Figure 1 are space (r and z in a cylindrical reference system) and time (t) dependent but we do not show this explicitly for the sake of a clear notation.Most features of the simulation framework have been detailed in [23]; nevertheless, here we summarize the main model ingredients and the extensions.
Intra-and extracellular fluids are described via the electrodiffusive Poisson-Nernst-Plank (PNP) transport model (see (1)(2)(3) in Figure 1) [45], implemented via a multiphysics coupling between the "electrostatics" and "transport of diluted species" interfaces from the AC/DC and electrochemistry modules [44].A reference electrode (RE in Figure 1) is placed to set the ground potential reference in the bulk of the extracellular fluid, where also the ion concentrations are at their baseline value.
The phospholipidic neuron membrane is described with a thin layer approximation according to (9) that saves a fine mesh.This enforces the continuity of the normal component of the displacement field vector across the intra-and extracellular boundaries and naturally creates diffuse layers at these interfaces.Ion channels, fundamental for generating action potentials (APs), are located at each mesh point of the lipidic layer.Their dynamics for AP genesis is described via the Hodgkin-Huxley (HH) model [9], [33] assuming a rest membrane potential Vr=-65 mV.The HH's gating variables are computed by solving a set of ordinary differential equations (ODE) at each membrane's mesh point (6).Consequently, the corresponding current density of the ion channels is updated through (7) [33] assuming Vr as resting potential.boundary flux to the electrolyte according to (10).The transmembrane stimulus to elicit action potentials is an additional Na + boundary flux across the upper portion of the neuron membrane (i.e., the one that does not face the electrode and the substrate), which is a simplified model for excitatory synapses [46].The membrane potential, V in (4), controls the non-linear conductances.
The reversal potentials for the ions are commonly included as constant parameters (i.e., set to their baseline values) in the ODEs.In our case instead, they can be dynamically updated by sampling the ion concentrations at each time step and assuming that changes in the concentrations instantaneously turn into changes in potential through (5).This last option entails that we self-consistently account for the impact of extracellular concentrations on the membrane potential.This topic is examined in subsection III-C.
In addition to our previous work [23], here we also account for uneven distributions of ion channels over the membrane.Indeed, different morphological regions of the same neuron (e.g., soma, dendrites, axon hillock, etc.) have intrinsically different channel densities [47]- [49] which might change locally if a neuron/electrode junction is present [34].As a support to this claim, fitting of extracellular recordings suggested accumulation/depletion of ion channels (e.g., Flowchart of the algorithm to derive the (b) elementary blocks of the lumped-elements equivalent circuit considering the dome-shaped neuron above to the planar sensing electrode in Figure 2.a.The circuit should be composed of at least one junctional (j), one side, and one upper (u) block.The junctional and lateral (l) blocks can be instantiated M and N times, respectively, depending on the desired partitioning degree of the system.It can be adapted to the cases in Figure 2.b and 2.c with minor adjustments that will be addressed in section III.
K + , Na + ) in rat hippocampal [50], cortical [30] and dorsal root [42] neurons, HEK293 cells [51], human myocardial cells [52].Furthermore, ion channels can drift-diffuse across the cell membrane [48] and different shapes of extracellular recorded signals are also observed across excitable cells without arborizations, such as chromaffin cells [53].To reproduce these putative uneven distributions of channels, we alter the maximum channel conductance g X i of ion species X i in (7) by multiplication for the scalar µ X i , similarly to [42].These parameters account for the accumulation (i.e., µ X i >1) or depletion (i.e., µ X i <1) of ionic channels at the membrane/electrode interface and the concurrent depletion (i.e., µ X i <1) or accumulation (i.e., µ X i >1) elsewhere in the membrane to maintain constant the overall number of ion channels (see (8)).
For the sake of computationally efficient calculations and differently from [23], we avoid meshing the bulk of the sensing electrode; instead, we collect the transduced electronic signal at the electrode surface (see (13) in Figure 1).The Stern layer on top of the electrode, included via a thin layer approximation, leads to the formation of a diffuse layer in the electrolyte (11) and to a displacement current entering the electrode (12).No red-ox reactions are implemented at the electrode surface.The electrode is connected to the input impedance of an ideal amplifier, although our approach can in principle include readout circuits of any complexity.
Figure 2 shows the structures considered in the FEM analysis, all featuring cylindrical symmetry.Physical properties and dimensions are respectively reported in Table I and Table II.

B. Methodology to derive a lumped-elements equivalent circuit
The lumped-elements equivalent circuit model is based on a few preliminary approximations.First, the electrolyte's reactance is assumed large compared to its parallel resistance, given the time scale of the action potentials.Second, intracellular compartments are assumed equipotential.Third, the electrodes are ideally polarizable.Fourth, we consider only axis-symmetric structures (i.e., the circular electrode is perfectly aligned with the neuron and are described in a cylindrical reference system).
Figure 3 illustrates the methodology to derive and combine the elementary blocks of the lumped-elements equivalent circuit.The procedure is described here for the dome neuron coupled to the planar electrode of Figure 2.a.It can be adapted to the cases in Figure 2.b and 2.c with minor adjustments that will be addressed in section III.
The first step is to partition the junctional portion of the neuron membrane (i.e., the one facing the electrode).This is represented by a ring with a small innermost radius of 1 nm and an outermost radius rse.We partition this ring in concentric subrings of equal pitch: where M is the number of the desired compartments for the junctional region.The use of rings of equal pitch is the one giving the best agreement between the FEM and the equivalent circuit, although we also considered rings with the same area still achieving a good mutual agreement (not shown).For each of the M rings, we identify the internal, r in,k , and external, r ex,k , radii.The area of the k-th ring is with k = {1, ..., M }, r in,1 =1 nm, and r ex,M =rse and we associate it with an HH circuit compartment.We found that 1 nm is a reasonable approximation for the innermost radius as long as r ex,1 ≫1 nm.Assuming to not consider the polarization effects of the proteinglycocalyx layer in the neuron/electrode cleft [54], each membrane compartment is in direct contact with the underlying portion of the electrolyte cleft (see step 2).In this latter case, careful observation of the FEM simulations shows the presence of a diffuse layer (DL) at the membrane's interface (i.e., ions pile up onto the membrane surface and their concentration decays according to the Debye-length [45]).One might thus expect a DL capacitance in between the HH circuit block and the sealing resistances representing the cleft.However, during neuron activity, the ion channels originate a mass transport path through the DL that effectively shortens the C DL at low frequencies.This effect translates into a DL resistor (R DL ) in parallel to the DL capacitance.We quantified R DL by means of a few transient simulations applying transmembrane ionic current steps (not shown) and we confirmed that the DL resistance (<1 nΩ•m 2 in our case) short-circuits any DL capacitance (whose reactance at 1 kHz is ≈184 µΩ•m 2 ).We will also show in Figure 4 that an equivalent circuit with DL capacitance in between the HH compartment and the cleft resistances generates sensed signals not consistent with the physically-based FEM results.For this reason, in the following, we will always assume that the parallel resistive path is highly conductive and consequently we remove the C DL at the membrane from the equivalent circuit.We point out that further polarization effects may stem from Stern layers at the neuron membrane, but there is no consensus on the necessity to account for them [32], [35] or not [36] in equivalent circuits.Moreover, there are no available strategies to implement them in FEM models while maintaining the ion concentrations close to the membrane within physiological ranges [23].Thus, we neglect their effects, in line with the other works [14].
The second step (step 2 in Figure 3) consists in evaluating the sealing resistance of the thin junctional electrolytic cleft between the electrode and the membrane of each of the M compartments.
Here we assume a radial current density propagation; consistently, for each compartment, we split the cleft's sealing resistance into two concentric contributions, each representing half of the in-plane area of the compartment.The first resistance, referred as internal, extends from r in,k to the radius r (A/2) k = (r 2 in,k + r 2 ex,k )/2 that divides into two equal parts the k-th ring's area and reads: The second resistive contribution, referred as external, goes from r (A/2) k to the end of the k-th ring, r ex,k , and is given by: R (ex) where σ el =1.43 S/m [35] is the electrolyte conductivity, k = {1, ..., M }, r in,1 =1 nm, and r ex,M =rse (consistently with ( 14)).We verified (not shown) that assigning R (in) j ,k and R (ex ) j ,k according to (15) and ( 16) rather than splitting the ring resistance into two parts at (r in,k + r ex,k )/2 lead to a much better agreement between the FEM model and the equivalent circuit.
The choice of ( 15) and ( 16), and the considerations on the membrane made above, allow us to connect the terminal of the kth HH compartment to the middle node of the underlying resistive net representing the k-th electrolytic cleft.This leads us to build the circuit in Figure 3.b by adding consecutive RC blocks with "T"-like topology.Notice that, only for the first block, the R (in) j,1 is missing since it does not carry any current (for symmetry reasons).
The third and fourth steps in Figure 3 account for the junctional and side portions of the electrode, respectively.By assuming an ideally polarizable metallic electrode (hence, inert in physiological fluid) [23], an electrical double-layer (Stern plus diffuse layers) builds up at the electrolyte-electrode interface, although more complex circuit models of the electrode surface can be found [55].This is modeled as an EDL capacitance per unit area, c EDL , computed as the series of c Stern and c DL , multiplied by the junctional electrode surface area of each k-th compartment in (14).The EDL capacitance is placed in between the resistances representing the electrolytic cleft, exactly below the corresponding HH compartment.The only exception is for the electrode sidewall (step 4 in Figure 3.a) where is connected to the sealing resistance representing the most peripheric part of the electrolytic cleft above the electrode.Notice that, differently from the C DL of the membrane, the electrode EDL capacitance should be retained, since there is no redox process nor mass transfer between electrode and fluid.
Step five in Figure 3 computes the area of non-junctional lateral membrane compartments (i.e., the ones just beside the junctional ones).Similarly to step 1, we partition this lateral ring in N concentric rings of equal pitch: where N is the number of the desired compartments for the lateral region.Thus, identifying N internal, r in,p , and N external, rex,p, radii, we compute the p-th HH compartment's area as: with p = {1, ..., N }, r in,1 =rse, and r ex,N =rneu.The sixth step determines the value of the internal and external portions of the sealing resistances that describe the electrolytic cleft of the lateral non-junctional portion of the neuron, following the same procedure as in step two (see (15)(16)).The only difference is that the cleft outside the electrode is thicker than the one above the electrode by a thickness tse (i.e., t cleft,l =t cleft,j +tse): At the last step, we represent the upper part of the neuron membrane (i.e., the one not facing the electrode or the planar substrate) as a unique compartment with half the surface area of an ellipsoid with semi-axis rneu and hneu, as reported in Figure 2.a, and neglecting rc for simplicity.Since this portion of the membrane faces the bulk electrolyte region connected to the reference electrode, the corresponding HH circuit is connected directly to the ground.

C. Validation of the equivalent circuit
In this section, the responses of the equivalent circuit to action potentials (APs) are compared to the FEM simulations of the structure in Figure 2.a with default parameters as in Table I and II.The APs are elicited by a transmembrane current pulse lasting 0.5 ms and with an amplitude of ≈0.22 nA applied to the upper HH compartment (i.e., 0.5 A/m 2 ).This choice avoids any overlap between the stimulus and the recorded extracellular signal [23], [35], [56].We account for the accumulation/depletion of channels in the junctional area (and the consequent depletion/accumulation of non-junctional channels) acting on the µ X i parameter for the sodium and potassium ions in (7) of Figure 1.We choose the same set of parameter values as in [42] (see Table II) and we sweep one parameter at a time w.r.t. the uniform channel distribution case (i.e., all µ X i =1).In this latter case we observe a very small sensed extracellular signal in the order of nanovolts, in agreement with experimental [29], [49] and simulation [34], [42] findings.This is because the capacitive and ionic currents essentially balance each other; thus, the fA-ranged net current that b)

V sens [μV]
-80 exits from adjacent HH compartments and flows through the cleft (whose total resistance is in the order of a few units of MΩ) is not capable to produce appreciable extracellular potential variations.4.a reproduce the typical AP waveform, and clearly show good agreement between the two modeling approaches.This waveform is the same for most of the case studies examined in this manuscript because any local alteration of the channel density next to the electrode (junctional membrane part) is compensated by an opposite alteration in the non-junctional part of the membrane.Therefore, it is not reported again in the following, except for the results in subsection III-C.Figure 4.b demonstrates that the circuit built following the steps in Figure 3 with M=4, N=4 matches very well the responses of the FEM model for all values of µ X i .The impact of the choice of M and N will be analyzed at the end of this section.All curves start with a pulse due to the transmembrane current stimulus.The pulse vanishes just before the onset of the extracellular signal activity at 0.5 ms.Different waveforms can be observed, and interpreted as follows: i) a depletion (black curves)/aggregation (red curves) of only the sodium channels at the sensing electrode makes the Na + current smaller/larger than the capacitive current in the fast rising phase of the AP.This non-negligible net current flows through the cleft and generates a voltage drop in the electrolyte surrounding the electrode.As a result, small positive/negative peaks can be observed around 0.8-0.9ms, respectively.ii) The depletion (blue curves)/aggregation (fuchsia curves) of the potassium channels or a complementary alteration (accumulation (red curves)/depletion (black curves)) of sodium channels at the sensing electrode has similar effects on the slow-falling phase of the AP, as visible by the negative/positive peaks at about 2 ms.This is consistent with the complementary roles of the Na + and K + channels.iii) As expected, close-to-zero (few nVs) voltage signals are attained if the channel distribution is uniform (green curves).Waveforms i-iii) are consistent with those reported in [42].iv) If the Na + and K + channels are simultaneously unbalanced in the same direction (µ K =µ Na =0.4 and 2, orange and violet curves), the early peak at 0.8 ms is further enhanced.These are two relevant conditions that resemble the typical biphasic shape of the positive/negative derivative of the intracellular AP as reported in experimental data [29].Also in this case, as in i-iii, the response of the circuit matches well that of the reference FEM model.
Figure 4.c compares the sensed signal when we include a DL capacitance in series to each HH block before the connection to the resistances of the cleft (i.e., R 3).As discussed in the previous section, the diffuse layer capacitance outside the membrane is shorted by a small resistance accounting for the transport of ions, an effect captured by the physical FEM model, not by the lumped-elements circuit used for this plot.Indeed, the inclusion of the DL capacitance alone (open symbols) causes large deviations from the FEM's results (filled symbols) and confirms that a short circuit should be inserted instead (see Figure 4.b), as will be done in the remaining of this work.II).To better spot the differences, we set µ Na =0.8 and µ K =1 which gives a large positive signal.The circuit with all the compartments better replicates the FEM response w.r.t. the case without the upper HH but this last compartment does not play an indispensable role in the qualitative description of the transient responses as the signal deviates by just 14% from the FEM curve for the chosen parameter set.
Figure 4.e investigates the effect of the number of HH compartments on the accuracy of the signal waveform.As in Figure 4.d, we consider µ Na =0.8, µ K =1.For M=8, N=8 the equivalent circuit response is at most 8% off to the FEM one, and converges to it for a larger number of compartments (not shown).Reducing M and N, the waveforms start to deviate from the reference FEM.Nevertheless, they maintain a good match down to the case M=2, N=1 (4 compartments total).This trend breaks down for the case with M=1, N=1.This is because the extracellular potential at the electrode/neuron junctional cleft should be distributed, by its nature, over the radial direction (given the similar sizes of the electrode and the neuron), while in this case it is lumped in a single compartment instead.This reduces the averaging effect of the electrode surface on the extracellular potentials.Thus, for comparable neuron/electrode sizes, a single-compartment approximation of the sensing part (i.e., a point contact model) overestimates the signals (at least by a factor of two in this case).

III. RESULTS
Having established the robustness of the equivalent circuit model over a wide range of parameter variations, in this section, we perform extensive parametric studies to show which are the most critical parameters affecting the shape of the recorded signals.Still, this analysis will be carried out with both the FEM model and the equivalent circuit one to further test the validity of the latter.

A. Planar Electrode
Figure 5 shows the simulated signals captured by a planar sensing electrode during an action potential.Different parametric studies are presented.In each study, all parameters that are not swept are kept at their default value reported in Table I and Table II.Unless otherwise stated, the equivalent circuit is built following Figure 3 with M=2, N=1 and considering two sets of ion channel depletion conditions to better visualize the curves: µ Na =0.8, µ K =1 and µ Na =1, µ K =0.8.
Figure 5.a reports the Vsens waveform for a cleft thickness ranging from 10 to 100 nm, and shows excellent agreement between the FEM and the equivalent circuit models in all cases.For thin clefts, the electrode is strongly sealed by the neuron membrane above, and the signal increases due to the large cleft resistance seen by the ionic currents through the cleft.For thick clefts, the signal is attenuated due to the large cleft conductance.
Figure 5.b investigates the optimum electrode size for maximum Vsens amplitude by sweeping the electrode radius from 2 to 15 µm (while the neuron radius is fixed at 10 µm).In all cases, we keep constant the area where the ion channel density is altered, A j , to the default value in Table I, regardless of the electrode size.We see that, as discussed in [23], a small electrode collects only a fraction of the extracellular signal activity due to its limited surface area and small sealing resistance.On the other hand, a large electrode (compared to the neuron size) seems even worse, since its uncovered portion, grounded by the electrolyte close to the RE, contributes to reduce the signal coming from the neuron.For the simulated geometry, the optimal electrode radius is about 0.5-1 times the neuron radius (i.e., rse=5-10 µm) that corresponds to the condition A j ≤A se,j ≤ πrneu 2 .
Figure 5.c examines the impact of different neuron sizes.The morphology remains a dome with constant hneu of 5 µm and rneu ranging from 6 to 50 µm.The smallest neuron delivers the smallest signal to the recording system.The largest positive/negative signals are recorded with 20/50 µm sized neurons.In the latter case, the first observed peak exceeds the scales of the graphs up to approximately 1 mV (not visible).Notice that in the rneu=50 µm case, the circuit has M=1 (hence it boils down to a point contact model) due to the (≈10 times) larger dimensions of the neuron compared to the sensing electrode.This result and the one in Figure 4.e give a useful indication of the minimum number of compartments required to build accurate equivalent circuits.
Figure 5.d investigates the effect of the input impedance of the readout focusing on two configurations taken from [36]: i) 100 GΩ resistor in parallel to 10 pF capacitor (i.e., mainly capacitive above ≈0.2Hz); and ii) 1 MΩ resistor in parallel to 9.45 nF capacitor (i.e., mainly capacitive above ≈17 Hz).These two configurations can be assumed as purely capacitive in the typical AP's spectrum range (≈1-1000 Hz) [5].Consequently, the sensed extracellular signal, which couples to the readout through the C EDLs at the electrode (see Figure 3), undergoes an almost pure capacitive divider with the readout capacitance.As a result, the sensed waveforms are undistorted and maximized in the small readout capacitance configuration owing to an advantageous capacitive divider.Notice indeed that the signals detected with the second configuration (red curves) are magnified 100 times for better visualization.
Figure 5.e reports the sensed voltage by sweeping the extension of the area with altered channel density.The signals are maximized when A j equals the top surface area of the sensing electrode, that is A j =A se,j =78.5 µm 2 in this case.Beyond this condition, we do not find any improvements, even if we alter the channel density in the whole bottom membrane (A j =314 µm 2 ). Figure 5.f shows the Vsens for the ellipsoidal neuron in Figure 2.b.This alternative neuron morphology is useful to analyze the effect of non-constant cleft thickness on the ability of our approach to build accurate equivalent circuits that replicate the FEM response even in complex geometries.We used a piece-wise constant thickness approximation for the clefts whose thickness increases as we approach the edge of the neuron.In this case, we set N=1; thus, we have a constant thickness for the non-junctional lateral compartment.Figure 5.f shows the recorded signal for all µ X i combinations.Also in this case the extracted equivalent circuit matches well the FEM results despite the simplification of curved clefts as boxes with constant thicknesses.

B. Mushroom Electrode
Mushroom electrodes (Figure 2.c) are a promising category of extracellular protruding electrodes for neural sensing [4], [5], [37].We report in this section a few parametric studies of this electrode morphology.
Figure 7 shows the responses of the FEM and equivalent circuit models for the structure sketched in Figure 2.c, where the cap of the mushroom is engulfed by the neuron.Figure 6 shows how the engulfment condition converts the planar disk cleft above a planar electrode into an equivalent (i.e., with the same area) semi-ellipsoid junctional cleft (i.e., the one in between the neuron and the mushroom cap).This curvature bends the current density streamlines from the center of the mushroom cap to its edge and might affect the signal sensed by the electrode.As a matter of fact, Figure 7.a shows the response if one approximates the cleft at the mushroom's cap as a planar disk (dashed lines) rather than a curved one (solid lines).This approximation, indeed, is inconsistent with the FEM.
To consider the curvature of the cleft when computing the resistances in the junctional blocks we followed the formula in [57] for concentric rings with a constant cleft at the mushroom's cap.In analogy with the planar electrode in steps 1-3 of Figure 3, we split the junctional neuron/electrode cleft surface into M rings of constant pitch on the planar disk cleft and equivalently on the semi-ellipsoidal cleft, as shown in Figure 6.Each k-th ring subtends a k-th angle in the mushroom, thereby we divide the π/2 rad bending angle into M subangles.In this procedure, the 2D profile of the semi-ellipsoidal cleft in Figure 6 is regarded as a semi-circumference with approximately the same perimeter; thus, a generic angle α w.r.t. the vertical direction (z-axis) be computed from the corresponding radius r on the planar disk according to Exploiting again the analogies with the planar electrode in Figure 3, every k-th sub-angle comprises an internal angle contribution from α in,k to α (A/2) k , and an external one from α (A/2) k to α ex,k .Doing so, the k-th internal and external junctional sealing resistances for  curved clefts can be computed as [57]: where tg(•) is the tangent function, α (A/2) k is the angle at which the area of k-th compartment divides into two, k={1,...,M}, α in,1 = 0 rad, and α ex,M = π/2 rad.The transient response of the equivalent circuit varying the cleft thickness and the electrode size is reported in Figure 7.b and 7.c, respectively.Both plots show an excellent agreement between the FEM results and the equivalent circuit description, demonstrating once more the robustness of the proposed partitioning of the junctional neuron/electrode area, for the determination of accurate lumpedelement model parameters.
In Figure 7.b we see that the effect of the cleft is similar for the mushroom and planar electrodes (see Figure 5.a) for the same reasons discussed in the previous section.
Figure 7.c analyzes the effect of the mushroom size.The cap radius spans from 1 to 4 µm for constant rneu=10 nm.In all cases, the neuron fully engulfs the cap at a constant cleft regardless of the electrode size; thus, the engulfed surface area corresponds to the external surface of the cap, denoted as A se,j .In this test, we consider that the junctional area with altered channel density: i) changes according to the actual engulfed area, i.e., A j =A se,j (solid lines); or ii) is constant at A j =5.2 µm 2 (dashed lines).The results suggest that a large mushroom electrode collects a larger useful signal than smaller ones only if the aggregation/depletion of channels spreads throughout the entire engulfed membrane portion.On the other hand, if the area where the channel density is altered remains constant, the signal is maximized as long as A se,j ≤A j .Similar behavior has been observed also for the planar electrode in Figure 5.b.

C. Effect of ion diffusion on the reversal potentials
As a final analysis, we use the proposed FEM model to compute, and self-consistently update at each time-step, the reversal potentials of the ions in the HH model according to the actual ionic concentrations at the HH membrane ((5) in Figure 1).We compare these results with the FEM simulations with reversal potentials fixed to     their baseline values (see Table I), both computed on the structure in Figure 2.a.We denote the former FEM implementation variable reversal potential (VRP) model and the latter constant reversal potential (CRP) model.
Figure 8 shows the AP responses computed with both approaches for cleft thicknesses of 50 nm and 300 nm, and µ Na =0.8, µ K =1 or µ Na =1, µ K =0.8.In particular, Figure 8.a shows that the intracellular potential waveform almost overlaps in the two models, except for a modest discrepancy during the hyperpolarization phase (after ≈ 3 ms) where the two curves are slightly spread apart.Figure 8.b shows the potassium concentration at the extracellular side of the bottom membrane (in correspondence of the symmetryaxis of the structure), [K + ]e (where e stands for extracellular), during the evolution of the AP in panel .a.The intracellular [K + ] is almost constant (not shown).Results at µ Na =0.8, µ K =1 and µ Na =1, µ K =0.8 exhibit very similar features.The strong confinement of the ions in the thin electrolytic cleft of 50 nm causes a larger increase of [K + ]e compared to the t cleft =300 nm case.The relative increase in [Cl -]e and [Na + ]e in the cleft w.r.t. the rest baseline (not shown) are less pronounced compared to one of [K + ]e and non-influential for our discussion.
Figure 8.c shows the time-varying potassium reversal potential obtained with (5) and K + concentration from panel .b.For 50 nm thick cleft, the V K + remarkably deviates from the rest value (≈ -77 mV [33]), which is not the case for the 300 nm thick cleft.The change in [K + ]e affects the sensed signal waveform in panels d-e.The [Cl -] and [Na + ] reversal potentials computed with (5) at the bottom membrane (not shown) remain almost constant at their baseline value (see Table I) during the AP. Figure 8.d and .edepict the sensed voltage signal for µ Na =0.8, µ K =1 and µ Na =1, µ K =0.8, respectively.We found that for thin clefts (50 nm) the sensed signals with the VRP significantly depart from the CRP case; whereas, at large clefts (300 nm), this difference disappears.This shift is due to the K + accumulation in the cleft and the consequent change of V K + from the baseline value.As a matter of fact, we found an almost perfect convergence of VRP and CRP models for t cleft ≥300 nm.Even at 100 nm the difference between the two approaches is quite small.The equivalent circuit results closely match the FEM with constant reversal potential since both models do not include corrections to the reversal potentials.This analysis suggests that including the sampling of the extracellular potassium concentration and self-consistent update of the reversal potential may be important when dealing with ultra-thin clefts.One should however consider that Figure 8 represents a kind of worst-case scenario, since we have not included other physiological mechanisms (e.g., Na + -K + -pumps) nor accounted for the presence of astrocytes surrounding the neuron that could buffer the extracellular potassium concentration.Furthermore, we assume that changes of [K + ]e instantaneously translate into changes of the reversal potential, although it is not obvious that this process may take place on the millisecond scale of an AP.These points definitely require further investigation in collaboration with the neuroscience community.

IV. CONCLUSION
We have developed a FEM simulation model for the analysis and optimization of extracellular recording electrodes that selfconsistently couples ion transport, action potential generation, and electrical signal collection.
The FEM framework is used to advance the state-of-the-art of area contact model by validating a robust and accurate approach to construct multi-compartment equivalent circuit models for the neuron/electrode system.Special care has been paid to the correct system partitioning and modeling of the connection between the HH blocks that describe the neuron membrane and the cleft resistance, and the rules to determine the minimum number of compartments for different neuron (dome and elliptical) and electrode (planar or mushroom) shapes.
A systematic analysis led us to derive or confirm design guidelines to maximize the quality of the sensed signal.The main findings are: 1) the electrolyte cleft between the neuron and the electrode should be as small as possible.For very thin clefts (below approximately 100 nm) ion diffusion inside the cleft affects the reversal potentials of the neuron membrane and perturbs the sensed signal; 2) the size of the electrode should be comparable to the size of the altered ion channels density portion of a neuron but should not exceed the neuron size; 3) the diffuse-layer capacitance at the membrane must not be included in the equivalent circuit since it is shorted by a very low resistive path created by the ion channels; 4) the minimum number of compartments needed for an accurate equivalent circuit description is 4 (M=2, N=1) for comparable neuron/electrode size, and 3 (M=N=1) for neurons at least ten times larger than the electrode.We highlight that, as commonly done in literature, this work assumed that the physical system has axial symmetry and thus the FEM simulations are solved in a cylindrical reference system to save mesh points.Additionally, this allowed us to adopt fairly simple analytical formulas for the lumped elements in the equivalent circuit based only on the radial and vertical dimensions of the physical system.Adaptations and extensions of the proposed method to domains without axial symmetry to account for misalignments between the electrode and the neuron are objectives of future work.
As a final remark, we point out that the developed FEM model is based on a general-purpose simulation framework; hence, it is upward scalable and amenable to include more complex descriptions of the neuron morphology and physiology, and more complete models of the readout circuit.These aspects, however, go well beyond the scope of this work.

Fig. 1 .
Fig. 1.Sketch of the modeling framework: neuron membrane, intracellular ("neuron") and extracellular ("electrolyte") fluids, sensing electrode, and readout circuitry.The boxes represent different physics, mutually coupled by the equations on the connecting arrows.All the physical variables are space (r, z in a cylindrical reference system) and time (t) dependent (not shown for improved clarity of the equations).The default values of the parameters and variables are given in Table I and II.
to M junctional-blocks up to N lateral-blocks 1 upper-block

μFig. 4 .
Fig. 4. Calculated intracellular action potential (V in ) (a) and corresponding sensed voltage (Vsens) (b-e) due to a transmembrane current stimulus for the default FEM model in Figure 2.a (filled symbols) and the corresponding equivalent circuit model (open symbols) built following the flowchart in Figure 3.(b) Signals sensed by the electrode for different µ X i .Good agreement between FEM and equivalent circuit is found with M=4, N=4, corresponding to 9 compartments in total.(c) Same as (b) but with the addition of C DL between the HH compartment and the sealing resistances.(d) Sensed signals for M=4, N=4, and µ Na =0.8, µ K =1 with and without the upper HH compartment.(e) Vsens computed for µ Na =0.8, µ K =1 and reducing the number of compartments from 17 (M=8, N=8) to 3 (M=1, N=1).In the last case, the circuit response strongly deviates from the FEM reference.Filled symbols=FEM model, open symbols=circuit model.

Figure 4
reports the intracellular (panel .a)and the sensed (panels .b-e)voltage signals for the default FEM model in Figure 2.a (filled symbols) and for the equivalent circuit model (open symbols).The intracellular signals in Figure

Figure 4 .
d depicts the signals at the sensing electrode considering or not the upper HH compartment in the equivalent circuit (values for junctional, non-junctional and upper areas are given in Table

2 Fig. 5 .Fig. 6 .
Fig. 5. (a-e) Simulated sensed signal Vsens for planar electrodes covered by the neuron as in Figure 2.a or (f) by the ellipsoidal-shaped neuron of Figure 2.b.The panels report the Vsens upon sweeping, once at a time: (a) cleft thickness; (b) electrode radius; (c) neuron radius; (d) input impedance of the readout; (e) area where channel density is altered; (f) neuron morphology.Filled symbols=FEM model, open symbols=circuit model, M=2, N=1.z

Fig. 7 .μ
Fig. 7. Simulated sensed signal Vsens for mushroom electrodes as in Figure 2.c using the FEM setup and the corresponding equivalent circuit (a) by approximating the mushroom's cleft as planar (dashed) and curved (solid), and, considering the curvature of the cleft in the circuit for different values (b) of cleft thickness and (c) of mushroom dimensions.Filled symbols=FEM model, open symbols=circuit model, M=2, N=1.

Fig. 8 .
Fig. 8. (a) Intracellular potential, (b) potassium ion concentration in the proximity of the bottom membrane, (c) corresponding potassium reversal potential, and (d-e) potential sensed AP transient response signals for the recording system in Figure 2.a.FEM simulations are performed with variable reversal potentials, i.e., (5) active (filled symbols) or constant reversal potentials (open symbols) for cleft thicknesses of 50 nm and 300 nm.The density of ion channels is: (b,c,d) µ Na =0.8, µ K =1 and (b,c,e) µ Na =1, µ K =0.8.
The main physical and geometrical parameters are reported in TableIand TableII, respectively.The simulation setup in (a) is used as the default case in the following.In the models, the electrode is connected by default to an ideal amplifier's input impedance Zamp=Ramp=100 GΩ.
Consistent with.b Set to to yield the same t cleft,l of the Figure 2.a.