How does the presence of neural probes affect extracellular potentials?

The cable equation [43–45] is of great importance in computational neuroscience, and it reads,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C_{\rm m}\frac{\partial v}{\partial t} + I_{{\rm ion}} =\eta\frac{\partial^2 v}{\partial x^2}, \label{Cable} \nonumber \end{align} \tag{ 12 }$$

where $v$ is the membrane potential of the neuron, $C_{\rm m}$ is the membrane capacitance, $I_{\rm ion}$ is the ion current density and $\newcommand{\e}{{\rm e}} \eta =\frac{h\sigma_{\rm i} }{4}$, where h is the diameter of the neuron, and $\sigma_{\rm i}$ denotes the intracellular conductivity of the neuron [43].

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 (12) by a discrete model [43]. 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,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_{\rm e}(x,y,z)=\frac{1}{4\pi \sigma_{\rm e} }\sum_{k}\frac{I_k}{\vert \mathbf{r}-\mathbf{r_k}\vert}. \label{CS} \nonumber \end{align} \tag{ 13 }$$

Here, $\sigma_{\rm e}$ is the constant extracellular conductivity (in all the implemented models, the milieu is assumed to be linear by using a constant $\sigma_{\rm e}$), $\mathbf{r_k}$ is the center of the kth compartment of the neuron, $\vert \mathbf{r}-\mathbf{r_k}\vert$ denotes the Euclidean distance from $\mathbf{r}=\mathbf{r}(x,y,z)$ to the point $\mathbf{r_k}$, and I_k denotes the transmembrane current of each compartment. This solution assumes that the extracellular milieu is purely conductive, infinite, and homogeneous. We denote this method as current summation approach (CS) [6].

As the silicon probes are made of insulated material, they could be approximated with the method of images (MoI) [12, 24, 25]. With the MoI the probe is assumed to be an infinite insulating plane, effectively increasing the extracellular potential by a factor of 2. 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. For the MoI, the summation formula (equation 13) reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_{\rm e}(x,y,z)=\frac{1}{2\pi \sigma_{\rm e} }\sum_{k}\frac{I_k}{\vert \mathbf{r}-\mathbf{r_k}\vert}. \label{MoI} \nonumber \end{align} \tag{ 14 }$$

As will be shown section 3.1, the peak scaling factor (1 and 2 for the CS and MoI solutions, respectively) of the modeled probes is modulated by the neuron-probe alignment, rotation, and by the probe type and it can be a value between 0 and 2 depending on these factors. Therefore, we also propose and compare a third current summation-based approach, namely scaled current summation (SCS), in which the scale factor is set to match the peak ratio between the hybrid solution (section 2.4.2) and the CS solution on the electrode with largest amplitude (e.g. 1.65 is used in section 3.3.1).

We implemented the same simulations presented in section 2.1 using the conventional modeling approach described above (CS) to compare them with the EMI simulations. We used LFPy [11, 12], running upon Neuron 7.5 [9, 10], to solve the cable equation and compute extracellular potentials using equation (13). As morphology, we used a ball-and-stick model with an axon with the same geometrical properties described in section 2.2. Similarly to the EMI simulations, we used a synaptic input in the middle of the dendritic region activated in the EMI simulation (z  =  360 $\mu$m) to induce a single spike and we observed the extracellular potentials on the recording sites. The synaptic weight was adjusted so that the extracellular largest peak was coincident in time with the one from the EMI simulation. 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⁻² so we set $g_L = g_{\rm pas}=0.06 \cdot 10^{-3}$ S cm⁻²) and we used an axial resistance R_a of 150 $\Omega$ cm⁻¹. The fixed_length method was used as discretization method with a fixed length of 1 $\mu$m, yielding 658 segments (23 somatic, 422 dendritic, and 213 axonal). Transmembrane currents were considered as current point sources in their contributions to the extracellular potential, following equation (13) (using LFPy pointsource argument of the RecExtElectrode class). The MoI and SCS solutions were calculated by multiplying the CS solution by a factor 2 and 1.65 (optimized scale factor using the hybrid solution).

The hybrid solution (HS) [21–23] combines the transmembrane currents for each neural segments computed with the cable equation and a finite element modeling for the extracellular space. The transmembrane currents are used as source terms in a finite element solution of the Poisson Equation in the extracellular space (equation (2), using an iterative solver for the Poisson problem, specifically, preconditioned conjugate gradients with algebraic multigrid preconditioning). With this approach, the probe can be explicitly modeled using insulating (Neumann) boundary conditions at the surface of the probe (equation (5)) and the differences between the HS and the EMI solution lie in differences regarding the modeling of the neuron dynamics, such as the self-ephaptic effect. The HS requires that a FEM simulation is run for each timestep of the transmembrane currents, each time setting the source terms with the currents at the specific timestep. This makes it computationally expensive, especially, for long simulations. Alternatively, one could run a single FEM simulation for each neural segment with a unitary test current and then use the potentials computed at the recording sites as a static map for summing the contribution of all currents at each timestep. The latter approach can be also computationally complex, as the number of segments in the multi-compartment simulation can be quite high and it would require to store in memory a large number of finite element solutions.

The hybrid solution is a good and widely used approach to model a non-homogeneous extracellular space, especially in the peripheral nervous system literature [21–23]. However, it requires to run a finite element simulation for every neuron simulation, as transmembrane currents are located in different positions for different neurons.

In order to overcome this issue, we designed the probe correction method (PC) that relies on the reciprocity principle [46] and the principle of superimposition (given the assumption of linearity of the milieu expressed in section 2.4.1). The reciprocity principle states that if a current I₁ in a position $(x_1, y_1, z_1)$ generates a potential u₁ in a second position $(x_2, y_2, z_2)$, then the same current I₁ placed in $(x_2, y_2, z_2)$ will result in a potential u₁ in $(x_1, y_1, z_1)$⁶. Using this principle, we first simulated with a finite element method the extracellular potential generated by a test current (1 nA) from each electrode i of a specific probe (e.g. Neuronexus) in any point of the extracellular space and define it as $u_i(x_i, y_i, z_i)$, where $(x_i, y_i, z_i)$ is the relative position with respect to the electrode i. Also in this case we used an iterative solver for the Poisson problem (preconditioned conjugate gradients with algebraic multigrid preconditioning). Then, leveraging on the reciprocity and superimposition principles, we mapped the contribution of each transmembrane current to the potential at each electrode i as: $u_{ik}=I_k u_i(x_k, y_k, z_k)$, where $(x_k, y_k, z_k)$ is now the relative position between the kth neural segment and the electrode i, and I_k is the transmembrane current for the kth neural segment. The potential at each electrode i can be computed as:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle u_i=\sum_ku_{ik}=\sum_k I_k u_i(x_k, y_k, z_k). \nonumber \end{align*}$$

The PC method allows to pre-compute the effect of a probe in the extracellular space and then use this mapping for any neural model, without the need to run a full FEM simulation. The number of FEM solutions that need to be computed and stored during the pre-mapping is equal to the number of electrodes in the probe.

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 3 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, while the orange traces show the potential when the probe is present in the extracellular space. In this case the probe tip was placed 40 $\mu$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 play an important role. In particular, for the Neuronexus probe the minimum peak without the probe is  −21.09 $\mu$V and with the probe it is  −41.26 $\mu$V: the difference is 20.17 $\mu$V. For the Neuropixels probe the peak with no probe is  −21.2 $\mu$V, with the probe it is  −44.36 $\mu$V and the difference is 23.16 $\mu$V. In case of the microwire type of probe, the effect is minimal: the minimum peak without the probe is  −16.85 $\mu$V, with the probe it is  −15.82 $\mu$V, and the difference is about 1.03 $\mu$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 $\mu$m), it extends for 30 $\mu$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.

Zoom In Zoom Out Reset image size

The MEAs, electrically speaking, are like insulating walls that do not allow currents to flow in. The insulating effect can be appreciated in figure 4, in which the extracellular potential at the time of the peak is computed in the [10, 100] $\mu$m interval in the x direction and in the [$-200, 200$] $\mu$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 4(C) shows that the effect at the tip of the MEA probe is almost null).

Zoom In Zoom Out Reset image size

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 $\mu$m to a maximum distance of 67.5 $\mu$m.

In figures 5(A)–(C) we plot the absolute peak values with (orange) and without probe (blue), as well as their 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 $\mu$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 $\mu$V (88.8% of the amplitude without probe) at 7.5 $\mu$m probe-membrane distance—and it decreases as the probe distance increases. At the farthest distance, where the probe tip is at 75 $\mu$m from the somatic membrane, the difference is 4.38 $\mu$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 $\mu$V (95.9% of the amplitude without probe) when the probe is closest and a minimum of 5.08 $\mu$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 $\mu$m from the somatic membrane. For the microwire, instead, 7.5 $\mu$m is the distance to the beginning of the cylindrical probe, whose tip extends in the x direction for 30 $\mu$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 (13)).

Zoom In Zoom Out Reset image size

In panel (D) of figure 5 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.

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 likely to be 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 $\mu$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 $\mu$m from the center of the neuron. In figure 6(A) 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 $\mu$m (the half width of the probe is 57 $\mu$m for Neuronexus and 35 $\mu$m for Neuropixels).

Zoom In Zoom Out Reset image size

In order to evaluate the effect of rotating the probes, we ran simulations with the probe at 70 $\mu$m distance (to accommodate for different rotations), coarse 2 resolution, boxsize 4, and rotations of 0, 30, 60, 90, 120, 150, and 180°. In figure 6(B) 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 4(C)), 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.

In order to compare the EMI simulations to conventional modeling, we built the same scenario shown in figure 2(B) (Neuronexus probe) using Neuron and LFPy, as described in section 2.4. As conventional modeling assumes an infinite and homogeneous medium, we compared the EAPs obtained by combining the cable equation solution (equation (12)) and the current summation formula (equation (13)) with the EMI simulations without the probe. The extracellular traces for the current summation approach (CS, red) and the EMI model (blue) are shown in figure 7(A). 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 $\mu$V, while the value for the CS is  −27.95 $\mu$V (the difference is 4.91 $\mu$V). This difference, which has been previously observed, is intrinsic to the EMI model [6], and can be due to self-ephaptic effects [6, 13–18]. 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.

Zoom In Zoom Out Reset image size

The hybrid solution (HS) uses currents computed with the cable equation and runs a FEM simulation of the extracellular space, including the probe. In figure 7(B) we show the extracellular potential of the EMI simulation with probe (orange) and the HS (green). Also in this case we observe that the EMI solution yields slightly smaller amplitudes with respect to the HS (EMI peak:  −42.6 $\mu$V; HS peak  −46.15 $\mu$V; difference: 3.55 $\mu$V) and these differences can be once again traced back to underlying differences of the neural solver.

After having shown that there are intrinsic differences between the EMI model and solutions based on the cable equation (CS, HS), we now compare two computationally less expensive strategies that could be used to account for the probe effect in modeling of extracellular potentials.

The MoI and SCS are attractive candidates due to their almost null computational cost, as they only multiply all values by a constant factor. The factor for infinite insulated planes, as described in section 2.4.1, is 2, but as shown in figures 5 and 6, for MEA probes it is somewhere between 0 and 2 depending on the neuron-probe lateral shift and rotation. In this scenario, the neuron is perfectly aligned with the probe and there is no rotation. The peak ratio for the SCS was computed by dividing the largest peaks of the HS and CS solutions and it was set to 1.65. In figure 8(A) the EAP from the HS (green), from the MoI (pink), and from the SCS with factor 1.65 (grey) are displayed. The MoI (pink) overshoots the estimation of the extracellular amplitudes (MoI peak:  −55.89 $\mu$V; HS  −46.15 $\mu$V; difference: 9.74 $\mu$V). The SCS solution, expectedly, results in the same amplitude as the HS on the electrode with the largest peak, as the scaling factor was computed using the actual peak ratio between the HS and the CS solution. However, there are some discrepancies between HS and SCS. Figure 8(B) shows the distribution of peak ratios of all the 32 electrodes with respect to the HS peaks. The CS, MoI, and SCS solutions display a range of values in the peak ratios, showing that the amplitude modulation of the electrodes is not a constant value. This can be traced back to the fact that a lateral shift of the neuron reduces the peak ratio (figure 6(A)): electrodes on the side of the probe yield a lower effect than the ones at the center of the probe. Due to this variability, a correction strategy based on a constant scaling will not be able to accommodate for this effect.

Zoom In Zoom Out Reset image size

The probe correction (PC) solution, based on the reciprocity principle (section 2.4.3), results in a solution perfectly coincident to the HS, at a much smaller computational cost (see table 5). In figure 8(B) the PC ratios are depicted as a vertical line at 1 because the peak amplitudes are exactly the same as the HS. The PC approach, in fact, pre-maps the effect of each electrode on the extracellular domain, effectively modeling in an efficient way the distribution of peak ratios observed when using the CS, MoI, and SCS methods.

In table 4 we summarize the comparison results, showing maximum, minimum, average peak ratios and the peak ratio distribution standard deviation for all the pairwise comparisons analyzed in this section.

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 [35], but it clearly cannot reach a level of detail of some more realistic morphologies, such as the reconstructed models made available by various initiatives [1–5]. We quantified the amplitude shift due to the probe in the extracellular domain ($\sim\!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 [49] 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, even if it has been shown that for normal electrodes' impedance the effect of conductive and equipotential recording sites is negligible [32].

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 min. The Neuron simulations described in section 2.4 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 $\Omega$ rather than $\Omega_e$ and $\Omega_i$ only, while the domain for the transmembrane potential $v$ is $\Gamma$, but the space for $v$ is setup on all facets of the mesh. For simplicity of implementation, the $v$ unknowns on facets outside of $\Gamma$ 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 5 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 5, 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.

The EMI framework, due to its computational requirements, is presently not an alternative to conventional modeling involving the cable equation (equation (12)) and the current summation formula (equation (13)). However, for specific applications, it can provide interesting insights. The hybrid solution combines 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 [21–23, 25, 47]. However, the HS is also computationally expensive and it increases in complexity with longer simulation durations. Similar considerations can be made if Boundary Element Methods (BEM) [50] are employed instead of FEM ones, even though they are less computationally intense then the current FEM formulation. One possible drawback of BEM solvers is that they could not accommodate for anisotropic conductivity, while FEM solvers could in principle solve meshes with non-homogeneous conductivity between surfaces [51].

Another much faster option could be using approaches based constant scaling, such as MoI and SCS. However, even correcting with a right factor smaller than 2, the these methods cannot account for the variability of peak ratios among the electrodes (figure 8(B)). Therefore, we suggested here the probe correction (PC) method, which combines a one-time finite element simulation to model how each electrode of a specific probe affects the extracellular domain, and then uses the reciprocity principle to compute the potential on the recording sites arising from transmembrane currents. We showed that this method is able to reach the HS accuracy at a much smaller computational time (table 5), which is also not strongly dependent on the simulation duration. Moreover, the time required to compute the probe specific mapping (PC (map)) and loading the FEM solutions in memory (PC (load)) could be further reduced by decreasing the mesh resolution. This possibility should be further investigated with a convergence analysis, similar to section 3.2 for the EMI model.