Improving Neural Simulations with the EMI Model

Mathematical modeling of neurons is an essential tool to investigate neuronal activity alongside with experimental approaches. However, the conventional modeling framework to simulate neuronal dynamics and extracellular potentials makes several assumptions that might need to be revisited for some applications. In this chapter we apply the EMI model to investigate the ephaptic eﬀect and the eﬀect of the extracellular probes on the measured potential. Finally, we introduce reduced EMI models, which provide a more computationally eﬃcient framework for simulating neurons with complex morphologies.


Introduction
In recent years, huge efforts and resources have been spent in computational modeling of neuronal activity. For example, the Blue Brain Project (18; 16)(https: //bbp.epfl.ch/nmc-portal/welcome) has constructed and distributed several hundreds of biophysically detailed cell models (multi-compartment models) from rat sensory cortex. A similar effort is being conducted by the Allen Institute of Brain Science, whose cell-type database (8) (https://celltypes.brain-map.org/) includes hundreds of cell models both from mice and even from humans. As the experimental data used become more comprehensive and available, these models are expected to become elaborated and more accurate in reproducing neuronal dynamics. However, the modeling framework which is commonly used to simulate these multi-compartment models makes several key assumptions that might be violated for certain applications.
The most widely used approach to simulate neuronal dynamics of neurons is the cable equation. The solution of this equation enables one to compute transmembrane currents for each of the model compartment. In order to simulate extracellular potentials, we can use volume conduction theory and sum the individual contributions of the currents to the electric potential at any point in space (6). Whereas the use of this modeling framework has been the gold standard to simulate neuronal activity for decades (19; 6), there are some important assumptions that need to be discussed: • A neuron is represented as a discrete set of nodes. Multi-compartment models split the neuronal morphologies into a discrete set of segments. Therefore, neurons are not represented as a continuum and this might affect the accuracy of the simulations. However, this assumption can be alleviated by using very small segments that can accurately replicate the neuronal complex geometry.
• Extracellular potentials are assumed to be constant. When solving the cable equation, the extracellular potentials outside the membrane are assumed to be constant. This assumption is harder to relax, as it prevents to include so-called ephaptic coupling in the simulation (9; 1). Ephaptic coupling refers to the effect of extracellular potentials on the neuronal dynamics. The use of the EMI model allows one to include these phenomena in the simulation, both to simulate the effect other neurons' activity or the same neuron's activity (self-ephaptic) has on the membrane potential.
• The extracellular space is assumed to be homogeneous. The most common approach to compute extracellular potentials arising from neuronal currents is to use volume conduction theory with the point-source or line-source approximations (6), which assume that the extracellular medium is homogeneous (in addition to linear, isotropic, and infinite). However, this assumption is clearly violated when using extracellular devices to record neuronal activity, which introduce a clear inhomogeneity in the extracellular space. Extracellular probes can be explicitly modeled in the extracellular space with the EMI model and they show to greatly contribute to the recorded signals (3).
In applications where any of the assumptions listed above may be violated the EMI model (10, (1.30)) provides a suitable modeling framework. In particular, the geometry of the neuron (and the extracellular space) is accurately represented. Here we will show how the model is convenient to study both the ephaptic coupling of neurons (Section 7.3) and the effects of extracellular probes on the recorded electric potentials (Section 7.4). However, the detailed representation of the geometries makes EMI much more computationally intense than the standard modeling framework (21; 3). This limits the complexity of the simulation mainly to simple neuronal morphologies, such as ball-and-stick neurons (3). In order to target more realistic morphologies Section 7.5 discusses the reduced EMI model where the (three-dimensional) intracellular space is lumped to a curve that is the centerline of the neuron.

EMI Simulations of Neurons using the neuronmi Python Package
Before discussing the applications let us first briefly introduce a Python package called neuronmi1, which has been used for the simulation results reported further. neuronmi provides a high-level application programming interface (API) to enable users to easily set up and run EMI simulations of neurons.
The workflow of the neuronmi package consists of two parts. First, a mesh needs to be generated. This is done with the generate_mesh function, that uses gmsh (7) as backend. With this function the user can choose different kinds of neurons to place in the mesh and optionally place a probe device in the extracellular space. Mesh resolution and sizes of the bounding box can also be adjusted, as well as parameters of the neurons and the probe. In the following code example, we create a mesh with a ball-and-stick neuron (bas) and a microwire probe. By default, the center of the soma is at (0,0,0) μm, the dendrite extends in the positive z direction and the axon in the negative z direction.
import neuronmi mesh_folder = neuronmi.generate_mesh(neurons='bas', probe='microwire') Once a mesh is generated, the EMI simulation can be invoked with the simulate_emi function, which implements the finite element method for the multi-dimensional mixed formulation (14, (5.14)) following the discretization proposed in (20). Through a set of parameters, the user can stimulate the neuron with a synaptic input, a step current, or a pulse current. Alternatively, the probe contacts can be used to stimulate the neuron extracellularly. By default, the neuron receives a synaptic input on its dendrite. The user can probe electric potentials u at any point in the mesh, while transmembrane currents i and membrane potentials v are available at facets on the neuron surface. The full solutions are also saved as pvd or xdmf files. The simulation is run as follows: Several parameters can be set to customize the mesh generation and the simulation. For further details, we refer to the package documentation (https://neuronmi. readthedocs.io/en/latest/).

Investigating the Ephaptic Effect between Neurons
Neurons are surrounded by the electrically conductive extracellular space. Groups of neurons create fluctuations in the local extracellular electrical field. These fluctuations in turn influence the intracellular electrical field through the ephaptic effect. Ephaptic coupling cannot influence neurons at rest, however, it can affect the spike timings of a neuron receiving suprathreshold stimulus.
We will illustrate how the EMI model can be used to compute the ephaptic coupling between two ball-and-stick neurons embedded in an extracellular space. The simulation is based on the neuronmi package detailed above. One of the neurons is stimulated with a synaptic input which elicits an action potential. The intracellular potential in the other neuron is sampled. We ran several experiments increasing the distance between the neurons. The intracellular potential is sampled in the centre of the ephaptically stimulated neuron (right) to measure the strength of the synaptic coupling. The deflection is 4.7 μV when the neurons are 5 μm apart and decreases to 3.2 μV when they are 40 μm apart with the soma of the stimulated neuron adjacent to the axon of the other neuron.
While in these simulations we only showed the effect of a single spike on an adjacent neuron, the occurrence of synchronous activity in populations of neurons can cause larger degrees of ephaptic coupling. The neuronmi package enables users to instantiate several neurons in the mesh and to define their inputs (which can be synchronous), hence allowing in principle to investigate ephaptic effects at the population level.

Investigating the Effect of Measuring Devices on Extracellular Potentials
While the presence of recording devices is usually ignored in the computation of extracellular potentials, recent findings suggest that newly developed silicon-based devices, or Multi-Electrode Arrays (MEAs), have a strong effect on the measured signals (3).
Using neuronmi, which is built on the EMI model, one can easily incorporate the neural devices in the mesh and investigate their effects on the recorded signals. To demonstrate this, we built meshes with a simple ball-and-stick neuron and different types of neural probes in its vicinity: Microwire: the first type of probe represents a microwire. For this kind of probes we used a cylindrical insulated model with 30 μm diameter. The extracellular potential, after the simulations, was measured at the center of tip of the cylinder. The microwire probe is shown in Figure 7.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 7.2B.

Neuropixels (MEA)
: the third type of probe model represents the Neuropixels silicon MEA (11). The original probe has more than 900 electrodes over a 1 cm shank. The probbe is 70 μm wide and 20 μm thick. In our mesh, shown in Figure 7.2C we used 24 12x12 μm recording sites arranged in the chessboard configuration with an inter-electrode-distance of 25 μm (11).
We ran EMI simulations using the meshes with and without the probe in the extracellular space (Figure 7. Microwire probes do not affect the recorded potentials, with an EAP peak of −21.63 μV without the probe and of −20.53 μV with the probe (Figure 7.2D). However, when recording with silicon MEAs, the extracellular potentials are strongly affected. For the Neuronexus probe (Figure 7.2D), the EAP peak without the probe is −30.47 μV, while with the probe it is −56.09 μV (peak ratio=1.84). For the Neuropixels probe (Figure 7.2E), the EAP peak without the probe is −32.73 μV, while with the probe it is −63.63 μV (peak ratio=1.94). The probe effect is probably due to the insulating properties of the silicon probes, which act as electrical walls for the generated currents. For further details and analysis we refer to (3).

Reduced EMI Model
In the examples considered thus far the problem geometry was simple allowing for computations on a serial desktop computer. In order to apply the EMI framework to realistic neurons two challenges need to be addressed: representation of neurons in the form of a finite element mesh and efficient solvers for large linear systems due to the EMI equations. However, even with the order optimal algorithms discussed in (12, Chapter 6) and efficient mesh generators for neuron surface geometries, see e.g. (17), the computational cost of the 3D-3D EMI models remains large (compared to the conventional approaches). As a computationally feasible alternative we shall next discuss the 3D-1D models.
Topological order reduction is a modeling technique used e.g. in reservoir simulations (4) or studies of tissue perfusion (5), which exploits geometrical properties of the system in order to derive its reduced model. Viewing a dendrite (branch) as generalized cylinder with length L and radius R we observe that R L and that in a typical domain of interest the neuron's volume is negligible compared to its surroundings. This property motivates a reduced representation of the neuron in terms of a (one dimensional) curve, the centerline, along with a function R, which provides radius of the crossection at each point of the line. An illustration of the concept can be seen in Figure 7.3. Thus, referring to the notation of Chapter 5, Ω i is reduced to a line while Ω e newly occupies the entire domain, i.e. Ω = Ω e . In turn, the reduced EMI model presents a coupling between three dimensional extracellular space and the one-dimensional intracellular space. We remark that the membrane is one-dimensional as well.
In order to apply order reduction to the EMI model we consider the singledimensional primal formulation (5.6). Note that therein, the coupling on the membrane Γ requires that both u e and u i are restricted from Ω e and Ω i respectively by the dedicated trace operator. In a reduced model Ω i = Γ, Γ = Λ and thus u i needs not to be restricted. On the other hand, restriction of u e to curve Λ can no longer be realized as a trace since such an operation is not well-defined for H 1 functions, see e.g. (5). To define a value of the extracellular potential on Λ let us introduce an averaging operator Here C R (x) = {y ∈ Γ; (y − x) ⊥ n(x)} with n(x) being the unit tangent vector of Λ at x, cf. Figure 7.3. Thus, Πu e is computed by sampling u e on the original (two- dimensional) neuron surface. However, in practical computations we assume that Γ has a circular cross section so that |C R (x)| = 2πR(x).
Using Π the reduced single-scale primal formulation of the EMI model reads: Find Here, the factors 2πR arise in reducing the integration domain from Γ to Λ and similarly for πR 2 and Ω i . Thus, defining the reduced specific membrane capacitance C m = 2πRC m and the reduced intracellular conductivity σ i = πR 2 σ i the operator form of (7.1) can be seen to be (6.2) with the new restriction operators T i = I and T e = Π. Note also that in (7.1) the function f , which characterizes the membrane dynamics, is defined on Γ.
For the proof of well-posedness of (7.1) as well as a detailed discussion of modeling assumptions, which allow for model reduction from (5.6) the reader is referred to (4). Moreover the reduced multi-dimensional primal formulation (5.8) is analyzed in (13).
To assess the reduced model (7.1) the surface mesh of a rat neocortex neuron RatS1-6-39 from the NeuroMorpho database (2) has been generated using AnaMorph (17). The neuron has been embedded into a (box-shaped) domain Ω such that |Ω|/|Ω i | ∼ 100 with the resulting geometry meshed by gmsh. The full 3D-3D single-dimensional primal formation has been used to compute the response to a 1 ms synaptic stimulus of 50 nA. Referring to Figure 7.4 the lower branch close to node number 4 has been stimulated. Using the centerline representation of the neuron the response has also been computed with (7.1). We remark that P 1 -P 1 elements were used with both formulations and that spaces were setup on conforming meshes. In particular, with (7.1) the discretization of Λ consisted of the edges of the cells of Ω. However, such a geometrically conforming discretization is not required in the reduced model. In fact, the meshes of Λ and Ω can be independent, see (4). The reduced model then resulted in 4587 unknowns, which is to be contrasted with 18248 unknowns due to (5.6). In turn, the simulation time using the reduced model is about 110 seconds while the full EMI model required cca. 340 seconds to complete.
The two models are compared in Figure 7.4 which shows values of the computed intracellular potentials at different points along the centerline. In general, there is qualitative agreement between the model predictions. However, the reduced model can be seen to tend to underestimate both the minima and the maxima, while the excitation occurs faster compared to the full model. In addition to the intracellular potentials, the extracellular potentials were compared by sampling 6.52μm away from the soma center (node 1 in Figure 7.4). We remark that the soma radius was 5.71μm. It can be seen that with −2.99μV < u e < 1.07μV for (7.1) and −1.99μV < u e < 0.67μV for (5.6) the reduced model overestimates the extrema. As with the extracellular potentials there is a temporal shift in the response; the negative peak is observed at 1.26ms, respectively 1.87ms.
While results of the comparison in Figure 7.4 shall be viewed as preliminary we argue that they illustrate sufficiently the potential of reduced EMI models. In particular, the reduced model is able to capture qualitatively the properties of the full EMI simulations. However, clear differences, especially in the temporal shifts of the peaks, have been observed. In the future we aim to investigate if suitable scaling of the stimulus and/or the parameters of the membrane ODEs can reduce the prediction error. In addition, the modeling error of the reduced model shall be evaluated similar to (15). More specifically, the soma, being approximated as a sphere in the 3D-3D model, cannot be represented as a slender cylinder (unlike the dendrites and axons). Thus the model reduction assumptions are not met on the soma. While localized in space, the effect of this error on temporal predictions should be analyzed. In turn, improved reduced models, which take into account the spherical nature of the soma, e.g. in construction of averaging operators, might be needed.

Conclusions
In this chapter, we have showcased some applications in which the EMI model can be a viable alternative to standard modeling frameworks in order to investigate aspects of neuronal activity and recorded signals.
We introduced an open-source software package named neuronmi to easily assemble meshes including simplified neurons and probes in the extracellular space. Finally, we performed preliminary investigations into accuracy and solution cost of a reduced 3D-1D EMI model suitable for simulating complex neuronal morphologies.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/ by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made. The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.