Can the presence of neural probes be neglected in computational modeling of extracellular potentials?

Objective: Mechanistic modeling of neurons is an essential component of computational neuroscience that enables scientists to simulate, explain, and explore neural activity. The conventional approach to simulation of extracellular neural recordings first computes transmembrane currents using the cable equation and then sums their contribution to model the extracellular potential. This two-step approach relies on the assumption that the extracellular space is an infinite and homogeneous conductive medium, while measurements are performed using neural probes. In general, large probes may be approximated as infinite non-conductive planes which will typically double the potential, but this estimate may be questioned for smaller probes. The main purpose of this paper is to assess to what extent the presence of the neural probes of varying shape and size impacts the extracellular field. Approach: We apply a detailed modeling framework allowing explicit representation of the neuron and the probe to study the effect of the probes and thereby estimate the effect of ignoring it. We use meshes with simplified neurons and different types of probe and compare the extracellular action potentials with and without the probe in the extracellular space. Main results: Our computations show that small probes (such as microwires) hardly influence the extra-cellular electric field and their effect can therefore typically be ignored. In contrast, Multi-Electrode Arrays (MEAs) significantly affect the extracellular field by magnifying the recorded potential. While MEAs behave similarly to infinite insulated planes, which double the potential, we find that the MEA gain factors strongly depend on the neuron-probe alignment and probe orientation. Significance: Ignoring the probe effect might be deleterious in some applications, such as neural localization and parametrization of neural models from extracellular recordings. Moreover, the presence of the probe can improve the interpretation of extracellular recordings, by providing a more accurate estimation of the extracellular potential generated by neuronal models.


Introduction
Huge efforts have been invested in computational modeling of neurophysiology over the last decades. This has led to the development and public distribution of a large array of realistic neuron models, for example from the Blue Brain Project (bbp.epfl.ch [1,2]), the Allen-Brain Institute brain cell database (celltypes. brain-map.org [3]), and the Neuromorpho database (neuromorpho.org [4,5]). As experimental data become available, these models become both more elaborate and more accurate. However, some of the assumptions underlying the most commonly used models may not allow the accuracy necessary to obtain good agreements between models and experiments. For instance, it was pointed out in [6] that assumptions underlying the classical cable equation and the associated method for computing the extracellular potential, lead to significant errors both in the membrane potential and the extracellular potential. In the present paper we investigate whether the classical modeling techniques used in computational neurophysiology are sufficiently accurate to reflect measurements obtained by different types of probes, such as microwires/tetrodes, and larger silicon Multi-Electrode Arrays (MEAs). Traditionally, these devices are not represented in the models describing the extracellular field, and our aim is to see if this omission introduces significant errors.
The most widely accepted and used modeling framework for computing the electrophysiology of neurons is the cable equation [7,8,9,10,11], which is used to find current and membrane potentials at different segments of a neuron. One straightforward and computationally convenient way to model the extracellular electric potential generated by neural activity is to sum the individual contributions of the transmembrane currents (computed for each segment) considering them as point current sources or line current sources [7,11] using volume conductor theory. Although this approach represents the gold standard in computational neuroscience, there are some essential assumptions that need to be discussed. First, i) the neuron is represented as a cable of discrete nodes and the continuous nature of its membrane is not preserved. Second, ii) when solving the cable equation, the extracellular potential is neglected, but the extracellular potential is computed a-posteriori. Third, and foremost, iii) when computing extracellular potentials, the tissue in which the neuron lies is modeled as an infinite medium with homogeneous properties. The validity of these assumptions must be addressed in light of the specific application under consideration. The first assumption i) can be justified by increasing the number of nodes in the model, but assumption ii) is harder to relax since it means that the model ignores ephaptic effects. Therefore, this assumption has gained considerable attention [12,13,14,15,16,17,6]. However, the main focus of the present paper is assumption iii). More specifically our aim is to study the effect of the physical presence of a neural probe on the extracellular signals. Can it be neglected in the mathematical model, or should it be included as a restriction on the extracellular domain? Specifically, is the conventional modeling framework, ignoring the effect of the probes, sufficient to yield reliable prediction of extracellular potentials?
In order to investigate this question, we have used the Extracellular-Membrane-Intracellular (EMI) model [18,6,19]. The EMI model allows for explicit representation of both the intracellular space of the neuron, the cell membrane and the extracellular space surrounding the neuron. Therefore, the geometry of neural probes can be represented accurately in the model. We have run finite element simulations of simplified pyramidal cells combined with different types of probes, such as microwires/tetrodes, and larger silicon Multi-Electrode Arrays (MEAs).
Our computations strongly indicate that the effect of the probe depends on several factors; small probes (microwires) have little effect on the extracellular potential, whereas larger devices (such as Multi-Electrode Arrays, MEAs) changes the extracellular potential quite dramatically, resembling the effect of a non-cundictive infinite plane in the proximity of the neuron. The effect, however, depends on the neuron-probe alignment and orientation. The results may aid in understanding experimental data recorded with MEAs [20], it may improve accuracy when extracellular potentials are used to parametrize membrane models as advocated in [21], and to localize neurons from MEA recordings [22,23].
The rest of the article is organized as follows: in Section 2 we describe the methods used throughout the paper, with particular focus on the EMI model ( §2.2), the meshes ( §2.3), the finite element framework ( §2.4), and the conventional simulations used for comparison ( §2.5). In Section 3 we present our findings related to the effect of probes of different geometry on the extracellular recordings ( §3.1), the variability of our simulations depending on geometrical parameters of the mesh ( §3.2), and we compare our findings with conventional modeling and experimental data ( §3.3). Finally, we discuss and contextualize the work in Section 4.

Methods
In this section we introduce the modeling used to investigate the effect of the probes on the extracellular potential. In particular we describe the classical cable equation and the associated current summation approach, and we describe the EMI-model and the associated finite element method to solve this model.

Conventional modeling: the cable equation and current summation approach
The cable equation [24,25,26] is of great importance in computational neuroscience, and it reads, Here, v is the membrane potential of the neuron, C m is the membrane capacitance, I ion is the ion current density and η = hσi /4, where h is the diameter of the neuron and σ i denotes the intracellular conductivity of the neuron [24]. This equation is used to compute the membrane potential of a neuron and the solution is commonly obtained by dividing the neuron into compartments and replacing the continuous model (1) by a discrete model [24]. In order to compute the associated extracellular potential, it is common to use the solution of the cable equation to compute the transmembrane currents densities in every compartment, and then invoke the classical summation formula, Here, σ e is the extracellular conductivity, r k is the center of the k−th compartment of the neuron, |r − r k | denotes the euclidean distance from r = r(x, y, z) to the point r k , and c k denotes the transmembrane current of each compartment (cf. [6]).

The Extracellular-Membrane-Intracellular model
The purpose of the present report is to estimate the effect of introducing a probe in the extracellular domain on the extracellular potential. This can be done using a model discussed in [27,28,18,29,6] referred to as the EMI-model. In the EMI-model the Extracellular space surrounding the neuron, the Membrane of the neuron and the Intracellular space of the neuron are all explicitly represented in the model. The model takes the form σ e ∇u e · n e = 0 at ∂Ω p , (6) n e · σ e ∇u e = −n i · σ i ∇u i def = I m at Γ, In the simplified geometry sketched in Figure 1, Ω denotes the total computational domain consisting of the extracellular domain Ω e and the intracellular domain Ω i , and the cell membrane is denoted by Γ. Note that u i and u e denote the intra-and extracellular potentials, and v = u i − u e denotes the membrane potential defined at the membrane Γ. The intra-and extracellular conductivities are given respectively by σ i and σ e and in this work we assume that the quantities are constant scalars. The cell membrane capacitance is given by C m , and the ion current density is given by I ion . I m is the total current current escaping through the membrane. The EMI model is here considered with grounding (Dirichlet) boundary conditions, i.e. u e = 0, on the boundary of the extracellular domain (∂Ω e ) while insulating (Neumann) boundary conditions, i.e. σ e ∇u e ·n e = 0 were prescribed at surface of the probe (∂Ω p ). Note that the latter is a suitable boundary condition also for the conducting surfaces of the probe [30,31]. The resting potential (see Table 1) is used as initial condition for v.

Meshes
In order to implement the EMI model described above, the computational domain was discretized by unstructured tetrahedral meshes generated by gmsh [32]. We used a simplified neuron model similar to a ball-and-stick model [33,34], with a spherical soma with 20 µm diameter -whose center is in the origin of the axis -an apical dendrite of length L d = 400 µm and diameter D d = 5 µm in the positive z direction and an axon of length L d = 200 µm and diameter D d = 2 µm in the negative z direction. Both the axon and the dendrites are connected to the soma via a tapering in the geometry. On the dendritic side, the diameter at the soma is 8 µm and it linearly reduces to 5 µm in a 20 µm portion. On the axonal side, the axon hillock has a diameter of 4 µm at the soma and it is tapered to 2 µm in 10 µm.
The neuron was placed in a box with and without neural probes to study the effect of the recording device on the simulated signals. We used three different types of probes: Microwire: the first type of probe represents a microwire type of probe (or tetrode). For this kind of probes we used a cylindrical insulated model with 30 µm diameter. The extracellular potential, after the simulations, was estimated as the average of the electric potential measured at the tip of the cylinder. The microwire probe is shown in Figure 2A alongside with the simplified neuron.
Neuronexus (MEA): the second type of probe model represents a commercially available silicon MEA (A1x32-Poly3-5mm-25s-177-CM32 probe from Neuronexus Technologies), which has 32 electrodes in three columns (the central column has 12 recording sites and first and third columns have 10) with hexagonal arrangement, a y-pitch of 18 µm, and a z-pitch of 22 µm. The electrode radius is 7.5 µm. This probe has a thickness of 15 µm and a maximum width of 114 µm, and it is shown in Figure 2B.
Neuropixels (MEA): the third type of probe model represents the Neuropixels silicon MEA [35]. The original probe has more than 900 electrodes over a 1 cm shank, it is 70 µm wide and 20 µm thick. In our mesh,  Figure 1: Sketch of the simplified neuron geometry and its surroundings. The intracellular domain is denoted by Ω i , the cell membrane is denoted by Γ, and the extracellular domain is denoted by Ω e . The boundary of the probe is denoted by ∂Ω p and the remaining boundary of the extracellular domain is denoted by ∂Ω e . The normal vector pointing out of Ω i is denoted by n i , and n e denotes the normal vector pointing out of Ω e . L and D are the length and diameter of neural segments, respectively, and D 1 is the diameter of the hillocks in correspondence of the soma. In our simulations, we consider three types of probe geometry (see Figure 2). Note that the probe interior is not part of the computational domain. shown in Figure 2C we used 24 12x12 µm recording sites arranged in the chessboard configuration with an inter-electrode-distance of 25 µm [35].
In order to evaluate the effect of the described probes depending on the relative distance to the neuron (x direction), we generated several meshes in which the distance between the contact sites and the center of the neuron was 17.5, 22.5, 27.5, 37.5, 47.5, and 77.5 µm. Note that these distances refer to the beginning of the microwire tip (which extends in the x direction for 30 µm and to the MEA y − z plane (for the MEA probes the recording sites do not extend in the x direction). When not specified, instead, the distance for the microwire probe was 25 µm, 32.5 µm for the Neuronexus MEA probe, and 30 µm for the Neuropixels probe (center of the probe tip at 40 µm).
To investigate if and how the bounding box size affects the simulation, since the electric potential is set to zero at its surface, we generated meshes with five different box sizes. Defining dx, dy, and dz as the distance between the extremity of the neuron and the box in the x, y, and z directions, the three box sizes were: Moreover, we evaluated the solution convergence depending on the resolution by generating meshes with four different resolutions. Defining r n , r p , and r ext as the resolutions/typical mesh element sizes for the neuron volume and membrane, for the probe, and for the bounding box surface, respectively, the four degrees of coarseness were: coarse 0: r n = 2 µm, r p = 5 µm, and r ext = 7.5 µm coarse 1: r n = 3 µm, r p = 6 µm, and r ext = 9 µm coarse 2: r n = 4 µm, r p = 8 µm, and r ext = 12 µm coarse 3: r n = 4 µm, r p = 10 µm, and r ext = 15 µm At the interface between two resolutions, the mesh size was determined as their minimum. Further, having instructed gmsh to not allow hanging nodes the mesh in the surroundings of the neuron and probe is gradually coarsened to r ext resolution.
For each of the mesh configuration with varying probe model, box size, and coarseness we simulated the extracellular signals with and without the probe in the extracellular space and sampled the electric potential at the recording site locations (even when the probe is absent).

Membrane model and finite element implementation
On the membrane of the soma and the axon, the ionic current density, I ion , is computed by the Hodgkin-Huxley model with standard parameters as given in [36]. On the membrane of the dendrite, we apply a passive membrane model with a synaptic input current of the form where g s (x) = g syn , for x in the synaptic input area, The parameters of the dendrite model are given in Table 1, and the synaptic input area is defined as a section of the dendrite of length 20 µm located 350 µm from the soma, as illustrated in Figure 1. The EMI model (3)-(9) is solved by the operator splitting scheme and the H(div) discretization proposed in [19]. In this scheme a single step of the EMI model consists of two sub-steps. First, assuming the current Parameter Value Parameter Value Table 1: Model parameters used in the simulations. The parameters of the Hodgkin-Huxley model are given in [36].
membrane potential v is known the ordinary differential equations (ODE) of the membrane model are solved yielding a new membrane state and the value of v. In the second step (9) discretized in time with I ion set to zero is solved together with equations (3)-(8) using the computed value of v as input. This step yields the new values of intra/extra-cellular potentials u i , u e and the transmembrane potential v. The H(div) approach then means that the EMI model is transformed by introducing unknown electrical fields σ i ∇u i and σ e ∇u e in addition to the potentials u i , u e and v. Thus more unknowns are involved, however, the formulation leads to more accurate solutions, cf. [19, section 3.]. In our implementation the ODE solver for the first step of the operator splitting scheme is implemented on top of the computational cardiac electrophysiology framework cbc.beat [37]. For the second step, the H(div) formulation of the EMI model, see [19, section 2.3.3], is discretized by the finite element method (FEM) using the FEniCS library [38]. More specifically, the electrical fields are discretized by the lowest order Raviart-Thomas elements [39] while the potentials use piecewise constant elements. The linear system due to implicit/backward-Euler temporal discretization in (9) and FEM is finally solved with the direct solver MUMPS [40] which is interfaced with FEniCS via the PETSc [41] linear algebra library.

Cable equation simulations
We implemented the same simulations presented above using the conventional modeling described in Section 2.1 to compare them with the EMI simulations. We used LFPy [11], running upon NEURON 7.4 [9,10], to solve the cable equation and compute extracellular potentials using Equation 2. As morphology, we used a ball-and-stick model with an axon with the same geometrical properties described in Section 2.3. Similarly to our simulations, we used a synaptic input in the dendritic region to induce a single spike and we observed the extracellular potentials on the recording sites. To model the spatial extent of the electrodes, we randomly drew 50 points within a recording site and we averaged the extracellular potential computed at these points [11]. We used the same parameters shown in Table 1 (note that in NEURON conductances are defined in S/cm 2 so we set g L = g pas = 0.06 · 10 −3 S/cm 2 ) and we used an axial resistance R a of 150 Ω/cm. The fixed_length method was used as discretization method with a fixed length of 1 µm, yielding 623 segments (21 somatic, 401 dendritic, and 201 axonal). Transmembrane currents were considered as current point sources in their contributions to the extracellular potential, following Equation 2 (using LFPy pointsource argument of the RecExtElectrode class).
We used the same approach to simulate a pyramidal cell model from the Neocortical Microcircuit (NMC) Portal of the Blue Brain Project [1,2] (https://bbp.epfl.ch/nmc-portal/welcome) in order to model the experiment performed by Neto et al. [20], in which pyramidal cells from layer 5 in rat cortex were recorded both with an extracellular MEA (Neuronexus A1x32-Poly3-5mm-25s-177-CM32 -same probe modeled in this work) and with a juxtacellular glass pipette. The relative position of the MEA and pipette was calibrated using a microscope and controlled using micro-manipulators so that the position of the pipette tip with respect to the extracellular probe was known.
We recreated in modeling the scenario in experiment 2014_11_25_Cell3.0 (form http://www. kampff-lab.org/validating-electrodes/) by placing a detailed model of a thick-tufted pyramidal cell from layer 5 (L5_TTPC1_cADpyr232_1) at the position of the pipette and with the alignment de- scribed in [20]: the Neuronexus MEA is tilted by -48.2 • along the y axis and the soma of the cell is located at 31 µm from the closest recording electrode. A 3D reconstruction of the neuron and the probe is shown in Figure 3A. We ran a 500 ms-long simulation, during which the neuron was stimulated with a constant clamp current to the somatic compartment which induced 10 spikes ( Figure 3B). We then averaged the extracellular spikes to obtain the mean waveform ( Figure 3C) and we compared the simulated EAP to the ones recorded experimentally. We focused on the amplitude of the action potential and we investigated if our findings on the effect of the probe could make predictions closer to the observed data.

Results
In this section we present results of numerical simulation which quantify the effect of introducing probes in the extracellular domain on the extracellular potential. We show how this effect depends on the distance between the neuron and the probe, their lateral alignment, and the probe rotation. Furthermore, we evaluate the numerical variability of the solutions, we compare with the conventional modeling scheme, and finally report CPU-efforts for the simulations.

The probe effect
3.1.1 The geometry of the probe affects the recorded signals The first question that we investigated is whether the probes have an effect and, if so, how substantial this effect is and if it depends on the probe geometry. In order to do so we analyzed the extracellular action potential (EAP) traces with and without placing the probe in the mesh. In Figure 4 we show the EAP with and without the microwire probe (A), the Neuronexus probe (B), and the Neuropixels probe (C). The blue traces are the extracellular potentials computed at the recording sites when the probe was removed (Without probe : u noprobe ), while the orange traces show the potential when the probe is present in the extracellular space (With probe : u wprobe ). In this case the probe tip was placed 40 µm from the soma center, we used a box of size 2 and coarse 2 resolution. It is clear that the probe effect is more prevalent for the MEA probes than for the microwire, suggesting that the physical size and geometry of the probe plays an important role. In particular, for the Neuronexus probe the minimum peak without the probe is −21.09 µV and with the probe it is −41.26 µV: the difference is 20.17 µV. For the Neuropixels probe the peak with no probe is −21.2 µV, with the probe it is −44.36 µV and the difference is 23.16 µV. In case of the microwire type of probe, the effect is minimal: the minimum peak without the probe is −16.85 µV, with the probe it is −15.82 µV, and the difference is about 1.03 µV (the peak without the probe is even larger than the one with the probe). Note that the values for the microwire are slightly lower than the MEAs because even if the microwire tip center is at the same distance (40 µm), it extends for 30 µm in the x direction, effectively lowering the recorded potential due to the fast decay of the extracellular potential with distance. The recording sites of the MEAs, instead, lie on the y − z plane, at a fixed distance.
The MEAs, electrically speaking, are like insulating walls that do not allow currents to flow in. The insulating effect can be appreciated in Figure 5, in which the extracellular potential at the time of the peak is computed in the [10, 100] µm interval in the x direction and in the [-200, 200] µm interval in the z direction (the origin is the center of the soma). Panel A shows the extracellular potential with the probe (Neuronexus) and panel B without the probe. The currents are deflected due to the presence of the probe, and this causes an increase (in absolute value) in the extracellular potential between the neuron and the probe, as shown in panel C, where the difference of the extracellular potential with and without probe is depicted. The substantial effect using the MEA probe probably also depends on the arrangement of the recording sites: while for the MEAs, the electrodes face the neuron (they lie on the y − z plane) and currents emitted by the membrane cannot flow in the x direction due to the presence of the probe, for the microwire, the electrode is at the tip of the probe (at z = 0, extending in the x − y plane -see Figure 2) and currents can naturally flow downwards in the x direction, yielding a little effect ( Figure 4C shows that the effect at the end of the MEA probe is almost null). The probe (white area) acts as an insulator, effectively increasing the extracellular potential (in absolute value) in the area between the neuron and the probe (panel C, blue colors close to the soma and red close to the dendrite) and decreasing it behind the probe of several µV. The effect is smaller at the tip of the probe (the green color represents a 0 µV difference).

The amplitude difference is constant with probe distance
In this section we analyze the trend of the probe-induced error depending on the vicinity of the probe. We swept the extracellular space from a closest distance between the probe and the somatic membrane of 7.5 µm to a maximum distance of 67.5 µm.
In Figure 6A-B-C we plot the absolute peak values with (orange) and without probe (blue), as well as their The difference is small even when the probe is close to the neuron. (B) Neuronexus MEA probe: maximum peak without probe (blue), with probe (orange), and their difference (green). The difference is large at short distances and it decays at larger distances. (C) Neuropixels MEA probe: maximum peak without probe (blue), with probe (orange), and their difference (green). Also for this probe the difference is large at short distances and it reduces at further away from the neuron. (D) Ratio between peak with and without probe for Neuronexus (red), the Neuropixels (blue) and the microwire probe (grey). The ratio is almost constant at different distances and the average value is 1.9 for the Neuronexus, 1.91 for the Neuropixels, and 1.05 for the microwire probe. difference (green) for the microwire (A), Neuronexus (B) and Neuropixels (C) probes. For the microwire (A), as observed in the previous section, the probe effect is small and the maximum difference is 1.97 µV, which is 10.1 % of the amplitude without probe, when the probe is closest. For the Neuronexus MEA probe (B), at short distances the difference between the peaks with and without probe is large -40.5 µV (88.8 % of the amplitude without probe) at 7.5 µm probe-membrane distance -and it decreases as the probe distance increases. At the farthest distance, where the probe is at 72.5 µm from the somatic membrane, the difference is 4.38 µV, which is 90.2 % of the amplitude without probe. For the Neuropixels MEA probe (C) the effect is in line with the Neuronexus probe, with a maximum difference of 41.07 µV (95.9 % of the amplitude without probe) when the probe is closest and a minimum of 5.08 µV, which is still 116.1 % of the amplitude without probe, when the probe is located at the maximum distance. Note that the peak amplitudes on the microwire probe are smaller than the one measured on the MEAs at a similar distances. At the closest distance, for example, the Neuronexus MEA electrodes lie on the y − z plane exactly at 7.5 µm from the somatic membrane. For the microwire, instead, 7.5 µm is the distance to the beginning of the cylindrical probe, whose tip extends in the x direction for 30 µm. The simulated electric potential is the average of the electric potential computed on the microwire tip and it results in a much lower amplitude due to the fast decay of the extracellular potential with distance (see Equation 2). In panel D of Figure 6 we show the ratio between the peak with probe and without probe depending on the probe distance for the Neuronexus (red), Neuropixels (blue), and the microwire (grey) probes. The ratio for the microwire probe varies around 1 (average=1.05), confirming that the probe effect can be neglected for microwire-like types of probe, due to their size and geometry. Instead, when a MEA probe is used, the average ratio is around 1.9 and its effect on the recordings cannot be neglected.

The probe effect is reduced when neuron and probe are not aligned
So far, we have shown results in which the neuron and the probe are perfectly aligned in the y direction, but the probe effect is probably affected by the neuron-probe alignment, since the area of the MEA probe (we focus here on the Neuronexus and Neuropixels MEA probes as the effect using the microwire is negligible) facing the neuron changes depending on the lateral shift in the y direction and probe rotation.
To quantify the trend of the probe effect depending on the y shift, we ran simulations moving the probes at different y locations (10, 20, 30, 40, 50, 60, 80, and 100 µm) and computed the ratios between the maximum peak with and without the MEA in the extracellular space. The simulations were run with coarse 2 resolution and boxsize 5 and the probe tip was at 40 µm from the center of the neuron. In Figure 7A we show the peak ratios depending on lateral y shifts. The ratio appears to decrease almost linearly with the shifts, from a value of around 1.8-1.9 when the probe is centered (note that the peak ratio slightly varies depending on resolution and size, as covered in Section 3.2) to a value of around 1.2 when the shift is 100 µm (the half width of the probe is 57 µm for Neuronexus and 35 µm for Neuropixels).
In order to evaluate the effect of rotating the probes, we ran simulations with the probe at 70 µm distance (to accommodate for different alignments), coarse 2 resolution, boxsize 4, and rotations of 0, 30, 60, 90, 120, 150, and 180 • . In Figure 7B the peak ratios depending on the rotation angle are shown. For small or no rotations (0, 30 • ) the value is around 1.7 (note that we always selected the electrode with the largest amplitude, which might not be the same electrode for all rotations). For a rotation of 90 • the peak ratio is around 1 (the probe exposes its thinnest side to the neuron) and for further rotations the probe's shadowing effect makes the peak with the probe smaller (as observed in Figure 5C), yielding peak ratio values below 1. These results demonstrate that the relative arrangement between the neuron and the probe play an important role in affecting the recorded signals.

Solution dependence on domain size and resolution
We generated meshes of four different resolutions and five different box sizes, as described in Section 2.3, in order to investigate how the resolution and the domain size affect the finite element solutions. Since we are mainly interested in how the probe affects the extracellular potential and we showed that only for MEA probes this effect is large, we focus on the extracellular potential at the recording site with the maximum negative peak. We used the Neuronexus MEA probe for this analysis and the distance of the tip of the probe was 40 µm (the recording sites plane is at 32.5 µm from the somatic center). The recording site which experienced the largest potential deflection was at position (32.5, 0, −13) µm, i.e. the closest to the neuron soma in the axon direction. For a deeper examination of convergence of the EMI model refer to [6]. For resolutions coarse 0 and coarse 1 the box sizes 4,5 and 5, respectively, were too large to be simulated.
In Table 2 we show the values of the minimum EAP peak with and without Neuronexus probe, their difference, and their ratio grouped by the domain (box) size and averaged over resolution. Despite some variability due to the numerical solution of the problem, there is a common trend in the peak values as the domain size increases: the minimum peaks tend to be larger in absolute values, both when the probe is in the extracellular space (from −40.12 µV for box size 1 to −43.09 µV for box size 5) and when it is not (from −20.64 µV for Box size V peak with MEA ( µV) V peak without MEA ( µV) Difference ( µV) Peak ratio  ), which forces the electric potential at the boundaries to be 0. For this reason, a smaller domain size causes a steeper reduction of the extracellular potential from the neuron to the bounding box, making the peak amplitude, in absolute terms, smaller. The peak difference with and without the MEA probe appears to be relatively constant, but the peak ratio tends to slightly decrease with increasing domain size for the same reason expressed before (from 1.95 for box size 1 to 1.82 for box size 5). The solutions appear to be converging for box sizes 4 and 5, but the relative error (difference between box 1 and box 5 values divided by the value of box 5) is moderate (6.89 % for the peak with probe, 12.95 % for the peak without probe, and 4.14 % for the peak ratio). Nevertheless, the peak ratio values obtained with larger, around 1.8-1.85, domain sizes should be a closer estimate of the true value. Table 3 displays the same values of Table 2, but with a fixed box size of 2 and varying resolution (Coarseness). The relative error (maximum difference across resolutions divided by the average values among resolutions) of the peak with the MEA is 3.3 %, without the probe it is 6.65 %, and for the peak ratio it is 3.53 %.
Because the main purpose of this work was to qualitatively investigate the effect of various probe designs and the effect of distance, alignment, and rotation on the measurements, we used resolution coarse 2 and box size 2, which represented an acceptable compromise between accuracy and simulation time. For investigating the effect of probe rotation and side shift we increased the box size to 4 and 5, respectively, to accommodate the position of the neural probe. Finally, in Section 3.3 we increased the resolution to coarse 0 and used box size 3 to obtain more accurate results for the comparison with the cable equation simulations.

Comparison with cable equation and experimental data
In order to compare the EMI simulations to conventional modeling, we built the same scenario shown in Figure 2B (Neuronexus probe) using NEURON and LFPy, as described in Section 2.5. As conventional modeling assumes an infinite and homogeneous medium, we compared the EAPs obtained combining the cable equation solution (Equation 1) and the current summation formula (Equation 2) with the EMI simulations without the probe. The distance between the neuron soma center and the probe tip was 40 µm (recording sites on the x = 32.5 µm plane), the box size was 3 and the resolution was coarse 0. The extracellular traces for the cable equation approach (red) and the EMI model (blue) are shown in Figure 8A. The EAPs almost overlap for every recording site, despite some differences in amplitude. On the electrode with the largest peak, the value for the EMI solution is −23.03 µV, while the value for the cable equation solution is −28.46 µV (the difference is 5.43 µV). This difference has been observed before, and is intrinsic to the EMI model [6]. Note also that the condition that forces the extracellular potential to zero at the boundary of the domain causes a steeper descent in the extracellular amplitudes, as discussed in Section 3.2. When compared to the EMI solution with the Neuronexus MEA probe in the extracellular medium ( Figure 8B), with a peak of −42.6 µV, the difference is expectedly larger, 14.13 µV, because the cable equation approach neglects the probe effect.
Finally, we tried to qualitatively explain some experimental results from Neto et al. [20] including the probe effect. Note that the distances reported in the experiment have an uncertainty of at least 10.5 µm (computed outside the brain tissue and hence not considering the mechanical mismatch between air and the brain tissue)

Coarseness
V peak with MEA ( µV) V peak without MEA ( µV) Difference ( µV) Peak ratio  Table 3: Solution variability depending on resolution (Coarseness). The columns contain the maximum peak with the Neuronexus (MEA) probe, without the probe, the difference and ratio of the amplitudes with and without probe. The values are computed with a box size 2.
[20] and this uncertainty might affect the interpretation of the data. Nevertheless, to our knowledge, the dataset at matter is the only available dataset that combines extracellular MEA recordings and co-registration of a micro-pipette using precise micro-motors. In experiment 2014_11_25_Cell3.0, the reported pipette-MEA distance is 31 µm and the EAP peak is −421 µV. Simulating a detailed pyramidal cell placed at the same distance (see Section 2.5 for details on cell model) yields an average peak of −92.47 µV, much lower than the recorded one. Even reducing the distance by the 10.5 µm uncertainty and accounting for a soma radius of about 7.5 µm, placing the soma center at 13 µm distance from the probe, yields an average peak of −334.82 µV, still lower than the recorded one of 421 µV. When we correct for the soma radius and we add an uncertainty of 5 µm so that the distance is 18.5 µm, we obtain a peak value of −211.86 µV. Including the probe effect and correcting this value by the 1.9 factor estimated in Section 3.1, we obtain a peak value of −402.53 µV, which is 95.6 % of the recorded value.
Although the uncertainty in measuring the position of the neuron makes this experimental validation qualitative, we believe that the probe effect could represent the missing link in explaining recordings from MEA probes.  Table 4: FEM system size, resolution (Coarseness), box size, mesh parameters (number of cells, number of facets, number of neuron facets, and vertices), and CPU time to solve the EMI model with no probe in the extracellular space for different resolutions (Coarseness) and domain sizes (Box size). Note that for coarse 2 and coarse 3 the resolution of the neuron (r n = 4 µm) is the same.

CPU requirements
Whereas the EMI formulation represents a powerful and more detailed computational framework for neurophysiology simulations, it is associated with a much larger computational load. The simulations were performed on an Intel(R) Xeon(R) CPU E5-2623 v4 @ 2.60GHz machine with 16 cores and 377 GB RAM running Ubuntu 16.04.3 LTS. Table 4 contains the coarseness, domain size, number of tetrahedral cells, number of mesh vertices, total number of triangular cells (facets), facets on the surface of the neuron, the system size for the FEM problem, and the time in second (CPU time) to compute the solution for meshes without the probe in the extracellular domain. We show the results without probes in the extracellular domain, as they are they are computationally more intense due to the fact that the volume inside the probe is not meshed (although the resolution on the probe surface is finer, the resulting system size without the probe is larger than with the probe). The CPU requirements and the time needed to run the simulation strongly depends on the resolution of the mesh: the problem with coarseness 3 and box size 3 takes around 1 hour and 20 minutes (system size=745 789), while for the same box size and coarseness 0, the time required is around 22 hours (system size=5 271 370). The domain size also strongly affects the mesh size and computation time. For example, for the coarse 2 resolution, with respect to box 1, box 2 is 1.83x slower, box 3 4.16x, box 4 8.33x, box 5 20.51x.

Discussion
In this article, we have used an detailed modeling framework -the Extracellular-Membrane-Intracellular (EMI) model [6,19] -to evaluate the effect of placing an extracellular recording device (neural probe) on the measured signals. We used meshes representing a simplified neuron and two different kind of probes: a microwire (a cylindrical probe with diameter of 30 µm) and Multi-Electrode Arrays (MEAs), modeling a Neuronexus commercially available silicon probe and the Neuropixels probe [35]. We quantified the probe effect by simulating the domain with and without the probe in the extracellular domain and we showed that the effect is substantial for the MEA probes ( Figure 4B-C), while it is negligible for microwires ( Figure 4A). The amplitude of the largest peak using the MEA probes is almost twice as large (∼ 1.9 times) compared to the case with no probe, and this factor is relatively independent of the probe distance ( Figure 6D), but it is reduced when the neuron and the probe are shifted laterally ( Figure 7A) or when the probe is rotated ( Figure 7B). Moreover, we discussed the effect of varying the mesh resolution and of the size of the computational domain. We also compared our finite element solutions to solutions obtained by solving the conventional cable equation, and found that the latter gave result very similar to the finite element solution when the probe was removed from the extracellular space ( Figure 8A). Finally, we analyzed a publicly available experimental dataset [20] and showed that the recorded EAP amplitude could not be predicted by conventional cable equation-current summation modeling, but including the probe effect resulted in more accurate predictions. Therefore, we suggest that the probe effect can be a key element in modeling experimental data obtained with MEA probes. However, clearly further analysis is needed to clarify this matter. At present the computational cost of the EMI model prevents simulations of neurons represented using realistic geometries.

Comparison with previous work
In this work we used a finite element approach [19] to simulate the dynamics of a simplified neuron and to compute extracellular potentials using the EMI model. The use of FEM modeling for neural simulations has been performed before [27,42,43,44,28], but mainly as an advanced tool to study neural dynamics and ephaptic effects. In Moffit et al. [42], the authors simulated, using the cable equation approach, a neuron at 65 µm from a shank microelectrode with a single recording site, and then used the currents in a finite element implementation of the extracellular domain, including the shank microelectrode. They found that the amplitude of the recorded potential with the shank was 77-100 % larger than the analytical solution, but the spike shape was similar to the analytical solution (Equation 2), in accordance with our results (Figures 8A-B). The effects using MEA probes and varying distances, lateral shifts, and probe rotations were not investigated. In Ness et al. [31], an analytical framework for in-vitro planar MEA using the method of images [45] was developed. A detailed neural model was simulated using the cable equation and transmembrane currents were used as forcing functions for a finite element simulation to validate the analytical solutions. In the in-vitro case, in which the MEA is assumed to be an infinite insulating plane, the authors showed that the insulating MEA layer affects the amplitudes of the recorded potentials, effectively increasing it by a maximum factor of 2, which can be analytically predicted by the method of images (MoI). Using the MoI, the factor 2 can be explained as follows: for each current source an image current source is introduced in the mirror position with respect to the insulating plane, effectively doubling the potential in proximity of the plane and canceling current densities normal to the plane.
In this study, we investigated how large the effect of commonly used in-vivo probes is using the advanced EMI modeling framework. Our results are in line with these previous findings and we also show that the geometry, in terms of size and alignment of the probe, plays a very important role. We show that large silicon probes can be almost regarded as insulated planes when the neuron is aligned to them (potential increased by factor ∼ 1.9) for large ranges of distances ( Figure 6D). An interesting effect following the reduction of the amplitude factor with lateral shifts ( Figure 7A) is that neurons not aligned with the probe will be recorded with a lower signal-to-noise ratio (SNR) due to the smaller amplitude increase, assuming that other sources of noise are invariant with respect to the probe location (such as electronic noise and biological noise from far neurons). This might bias neural recordings towards identifying neurons that are closer to the center of the probe, rather than the ones lying at the probes' sides. However, this conclusion is speculative and might be affected by other factors, such as the distribution of neurons around the probe and their morphology (which contributes to the EAP). Therefore, ground truth information about the position of the recorded neurons and their reconstructed morphologies are needed for a quantitative evaluation of this phenomenon.

Limitations and extensions 4.2.1 Mesh improvements
The EMI model is, in principle, able to accurately represent the neuron and the neural probe. However, the accuracy of the model comes at the cost of computational resources. In order to be able to run simulations in a reasonable amount of time, the geometry of the neuron needed to be simplified considerably. First, we used a simple neuron in terms of a ball-and-stick with axon. This model is able to describe certain aspects of the neuronal dynamic [34], but it clearly cannot reach a level of detail of some more realistic morphologies, such as the cell model used in Section 2.5 and shown in Figure 3A. We quantified the amplitude shift due to the probe in the extracellular domain (∼ 1.9 on average for the MEA probes when neuron and probe are aligned), but this factor most likely also depend on the specific cell morphology that we used, and not only on the probe design and geometry. Therefore, we aim at extending the framework [46] for generating finite element meshes from publicly available realistic morphologies [5], allowing us to explore the probe effect for more complex morphologies.
Furthermore, we assumed ideal recording sites with an infinite input impedance which does not allow any current to flow in. In reality, recording electrodes have a high, but not infinite impedance that could be modeled by considering electrodes as an additional domain with very low conductance [30].

Computational costs
In Section 3.4 we showed that the EMI model is much more computationally demanding than conventional modeling using cable and volume conduction theory. For the simplest simulation performed in this study (coarse 3 and box size 1), a system with 337 515 unknowns was solved in about 40 minutes. The NEURON simulations described in Section 2.5 took ∼ 0.59 s to run, about 2400 times faster than the simplest EMI simulation performed here. However, because of our implementation and solution strategy for FEM , this factor should be considered as a rather pessimistic upper bound. In particular, the employed version of FEniCS (2017.2.0) does not allow for finite element spaces with components discretized on meshes with different topology. For example, the extra/intra-cellular potentials are defined on the entire Ω rather than Ω e and Ω i only, while the domain for the transmembrane potential v is Γ, but the space for v is setup on all facets of the mesh. For simplicity of implementation, the v unknowns on facets outside of Γ are forced to be zero by additional constraints and are not removed from the linear system. The LU solver thus solves also for the unphysical/extra unknowns and the memory footprint and solution times are naturally higher. The number of unphysical unknowns can be seen in Table 4 as a difference between total number of facets in the mesh and the number of facets on the surface of the neuron. For example, in the largest system considered here, avoiding the unphysical unknowns would reduce the system size by about 2 million.
In addition to assembling the linear system with only the physical unknowns, a potential speed up could be achieved by employing iterative solvers with suitable preconditioners. That is, fast PDE solvers for diffusion equations typically use around 1s per million degrees of freedom. As we here employ a H(div) formulation, we expect the solution to be computed in around 5 s per million degrees with multilevel methods. As shown in Table 4, 500 timesteps of solving systems with around one million degrees of freedom takes 82 600 s, which means 165 s per time steps. Hence, we may expect to speed up the solving procedure by around a factor 30 with better solvers. If further speed-up is required then finite element based reduced basis function method provides an attractive approach that should be addressed in future research.

EMI is not alternative to the conventional cable equation
The EMI framework, due to its computational requirements, is presently not an alternative to conventional modeling involving the cable equation (Equation 1) and the current summation formula (Equation 2). However, for specific applications, it can provide interesting insights. A hybrid solution could combine the cable equation solution to finite element modeling, in practice solving the FEM problem only for the extracellular space and using the transmembrane currents computed by the cable equation as forcing functions [42,31]. Similarly to what was suggested in Section 3.1, another much faster option could include a previous estimation of a probe-specific correction factor to accommodate for the amplitude mismatch. As a large MEA almost behaves like an insulating plane, also the method of images could straightforwardly be applied to correct for potential amplitudes (giving rise to a factor 2 [45,31]). While this is an interesting and simple option, it raises issues regarding the estimation of such correction factors, especially given the variability depending on relative lateral shift between the neuron and the probe ( Figure 7A) and probe alignment ( Figure 7B).

Significance of the probe effect
The effect of the recording device has not been fully taken into consideration in mathematical models of the extracellular field surrounding neurons. The probe effect needs to be considered when modeling silicon MEA, whose sizes are significantly larger than the recorded neurons. The assumption of an infinite and homogeneous medium is in fact largely violated when such bulky probes are in the extracellular space in the proximity of the cells. Although the tissue can be regarded as purely conductive and with a constant conductivity [47], these probes represent clear discontinuities in the extracellular conductivity, which strongly affect the measured potential due to their insulating properties. While the probe effect is large for MEAs, we found that it was negligible for microwire-type of probes, mainly for two reasons: first, the microwire is thinner and overall smaller than the MEA; second, the electric potential is sampled at the tip of the probe and in the entire semispace below the microwire currents are free to flow without any obstacle. When dealing with silicon MEAs, though, this effect could be crucial for certain applications that require to realistically describe recordings. For example, Gold et al. [21] used, in simulation, extracellular action potentials (EAP) to constrain conductances of neuronal models. Clearly, neglecting the probe effect would result in an incorrect parametrization of the models in this case.
Another example in which including this effect could be beneficial is when EAP are used to localize the somata position with respect to the probe. This is traditionally done by solving the inverse problem: a simple model, such as a monopolar current source [48,49,50], a dipolar-current source [48,51,52], line-source models [53,54], or a ball-and-stick model [55], is moved around the extracellular space to minimize the error between the recorded potential and the one generated by the model. Ignoring the probe might result in larger localization errors.
Recently, we used simulated EAP on MEA as ground truth data, from which features were extracted to train machine-learning methods to localize neurons [22,23] and recognize their cell type from EAPs [23]. When training such machine-learning models on simulated data and applying them to experimental data, neglecting the probe effect could confound the trained model and yield prediction errors.
Moreover, as we qualitatively showed in Section 3.3, explaining experimental recordings on MEA without considering the probe, like the data presented in [20], might cause discrepancies between the modeling and experimental results hard to reconcile. On the other hand, in order to fully explain and validate our findings, an experiment with accurate co-location of extracellular recordings and cell position (and ideally morphology) is required. For example, an experimental setup in which planar a MEA is combined with two-photon calcium imaging [56] could provide an accurate estimate of the relative position between the neurons and the MEA.
In conclusion, we presented numerical evidence that suggests that the probe effect, especially when using Multi-Electrode silicon probes, affects the way we model extracellular neural activity and interpret experimental data and cannot be neglected for specific applications.