Recorded potentials in the extracellular space (ECS) of the brain is a standard measure of population activity in neural tissue. Computational models that simulate the relationship between the ECS potential and its underlying neurophysiological processes are commonly used in the interpretation of such measurements. Standard methods, such as volume-conductor theory and current-source density theory, assume that diffusion has a negligible effect on the ECS potential, at least in the range of frequencies picked up by most recording systems. This assumption remains to be verified. We here present a hybrid simulation framework that accounts for diffusive effects on the ECS potential. The framework uses (1) the NEURON simulator to compute the activity and ionic output currents from multicompartmental neuron models, and (2) the electrodiffusive Kirchhoff-Nernst-Planck framework to simulate the resulting dynamics of the potential and ion concentrations in the ECS, accounting for the effect of electrical migration as well as diffusion. Using this framework, we explore the effect that ECS diffusion has on the electrical potential surrounding a small population of 10 pyramidal neurons. The neural model was tuned so that simulations over ∼100 seconds of biological time led to shifts in ECS concentrations by a few millimolars, similar to what has been seen in experiments. By comparing simulations where ECS diffusion was absent with simulations where ECS diffusion was included, we made the following key findings: (i) ECS diffusion shifted the local potential by up to ∼0. 2 mV. (ii) The power spectral density (PSD) of the diffusion-evoked potential shifts followed a 1/f2 power law. (iii) Diffusion effects dominated the PSD of the ECS potential for frequencies up to several hertz. In scenarios with large, but physiologically realistic ECS concentration gradients, diffusion was thus found to affect the ECS potential well within the frequency range picked up in experimental recordings.
The number of ions exchanged between neurons and the extracellular space (ECS) during a brief period of activity (i. e. , due to the integration of synaptic input and generation of a few action potentials) is typically too small to evoke significant changes in extracellular ion concentrations. In models of short-term electrical signalling of neurons, the ion concentrations of the main charge carriers (e. g. , K+, Na+, Cl-) are therefore commonly assumed to remain effectively constant. This assumption often holds also at longer time scales, due to the work done by neuronal and glial uptake mechanisms in maintaining ion concentrations close to baseline levels. However, during periods of intense neural signalling, the uptake mechanisms may fail to keep up, and ion concentrations in the ECS may change by several millimolars [1–5]. For example, the extracellular K+ concentration can increase from a typical baseline level of around 3 mM and up to levels between 8 and 12 mM during non-pathological conditions [4,6–8]. Ion-concentration shifts in the ECS will change neuronal reversal potentials and firing patterns [9–12], and too large deviations from baseline levels can lead to pathological conditions such as hypoxia, anoxia, ischemia, epilepsy and spreading depression [9,13–15]. One of the most common experimental methods for investigating neural activity is the measurement of electrical potentials with extracellular electrodes. Commonly, it is assumed that extracellular potentials predominantly reflect transmembrane cellular current sources, including synaptic currents and currents through active and passive membrane mechanisms in neurons and glial cells [16–19]. However, in scenarios where ECS concentration gradients become sufficiently large, electrical currents carried by diffusing ions in the ECS could in principle also give measurable effects on the extracellular electrical potentials (cf. , liquid junction potentials [20–22]). In support of this, local ion-concentration changes in the ECS are indeed often accompanied by slow local negative potential shifts, which can be on the order of a few millivolts [1,3, 13,23–27]. Whereas K+ buffering currents through the glia-cell membranes are believed to be the main source of these slow potential shifts [3,7], it has been estimated that also diffusive currents along extracellular concentration gradients could contribute by shifting ECS potentials by up to 0. 4 mV [3]. As ion concentrations in the ECS typically vary on the time scale of seconds [3,4, 28], it is nevertheless a priori unclear whether diffusion-evoked potential shifts would be picked up by the electrode measurement systems applied in most experiments, which typically have cut-off frequencies of about 0. 1–0. 2 Hz or higher (see e. g. , [29,30]). In most computational studies of ECS potentials, diffusive currents in the ECS are assumed to be negligible compared to the currents propelled by the electrical field (hereby termed field currents). This is, for example, an underlying assumption in volume-conductor theory which has been the basis for estimating ECS potentials from cellular current sources [18,31–36], and in estimation of current-source density (CSD) which predicts transmembrane neural current sources from recordings of extracellular potentials [29,34,37–40]. Another series of theoretical studies have aimed to incorporate possible effects of diffusion in the complex impedance environment of the extracellular medium [41–44], and have suggested that such effects may account for the 1/f-scaling observed for the LFP-power spectrum at low frequencies [41]. In neither of the above mentioned studies, however, ionic diffusion was explicitly modelled. The reason why diffusive effects are often neglected in models of extracellular fields, may be that the task of modelling it is challenging. This is because the study of diffusion requires an explicit tracking of all present ions and their spatiotemporal dynamics: i. e. , keeping track of only the electric currents and net electric charges is not sufficient. Existing electrodiffusive models have typically been based on the Poisson-Nernst-Planck (PNP) formalism [45–51]. The PNP formalism explicitly models charge-relaxation processes, which occur at spatiotemporal scales on the order of nanometers and nanoseconds. This requires an extremely high spatiotemporal resolution, which makes PNP models computationally expensive and unsuited for predictions at the tissue/population level [52]. However, a series of modelling schemes have been developed that circumvent the charge relaxation processes, essentially by replacing Poisson’s equation by the constraint that the bulk solution is electroneutral [28,52–58]. The electroneutrality condition is a physical constraint valid at a larger spatiotemporal scale, and thus allows for a dramatic increase in the spatial and temporal grid sizes in the numerical simulations. One of these simpler models were previously developed by our group [28,57], and is here referred to as the Kirchhoff-Nernst-Planck (KNP) scheme. The KNP scheme is a means of deriving the local potential in the intra- and extracellular bulk solution from the constraint that Kirchhoff’s current law should be fulfilled for all finite volumes (the sum of currents into a finite subvolume of bulk solution should be zero). In the current work we have developed a hybrid modelling formalism that allows us to compute electrodiffusive ion dynamics in the ECS surrounding active neurons. The formalism is briefly summarized in Fig 1. First, it utilizes the NEURON simulator [59,60], which is a standard tool for simulating morphologically complex neurons, to simulate the activity of a neural population and its exchange of ions with the ECS (Fig 1A). Second, it utilizes the KNP formalism [28,57] to compute the dynamics of ion concentrations and the electrical potential in the ECS surrounding the neurons (Fig 1B). The KNP scheme accounts for all electrical currents entering an ECS subvolume in the system (i. e. , transmembrane ionic currents, transmembrane capacitive currents, diffusive currents through the ECS, and field currents through the ECS), as well as for concentration-dependent variations in the ECS conductivity (see Methods). It computes the ECS potential from the constraint that all currents into a ECS subvolume should sum to zero (Fig 1C). In this way, the KNP-scheme accounts for effects of ionic diffusion on the ECS potential, and thus differs from previous simulation schemes for computing ECS potentials based on output from standard neuron simulators such as NEURON (e. g. [61]). We have here used the hybrid scheme to model a small tissue element consisting of a population of ten pyramidal neurons embedded in ECS (Fig 1A). Motivated by the layered structures of cortex and hippocampus, we assumed lateral homogeneity, so that all spatial variation occurred in the vertical direction. As neuronal model, we used a well established multicompartmental model of pyramidal cells [62]. The ten neurons were centered at the same depth level of the tissue. Although the choice of neuronal model was somewhat arbitrary, and the small tissue element was too simple to represent any particular biological system, the model gave rise to biologically realistic variations in ECS concentrations, and we regard it as a meaningful scenario for which we could explore how diffusive currents in the ECS can influence extracellular potentials. In simulations that evoked large, but not pathologically large, concentration gradients in the ECS, we found that diffusion gave rise to a detectable 1/f2 power law in the low-frequency part of the power spectral density (PSD) of the ECS potential. Furthermore, we found that diffusion influenced the PSD for frequencies as high as 1–10 Hz. This quantitative prediction was, of course, specific to the particular model setup used here. Although the relative effects of diffusion may be smaller in many realistic, more complex scenarios (see Discussion), we regard our findings as an important demonstration that in general, diffusive currents can not by default be assumed to have a negligible impact on ECS potentials. The article is organized as follows: In the Results section, we use the KNP scheme to explore the role of diffusive currents on electrical potentials in the ECS surrounding a population of pyramidal neurons. In the Discussion section, we discuss possible implications that our findings will have for the interpretation of data from extracellular recordings. The Discussion also includes an overview of the assumptions made in the presented model, and on how the framework can be expanded to allow for more thorough investigations of concentration-dependent effects on ion dynamics in neural tissue. A detailed derivation of the KNP-formalism is postponed to the Methods section (which is found at the end of the article).
Ten pyramidal neurons were simulated by running ten independent simulations on a single neuron model. As neuron model, we used a well established model developed for cortical layer 5 pyramidal cells [62]. Each neuron was driven by uncorrelated Poissonian input spike trains (with the same statistics for all neurons) through 10,000 synapses. Synapses were uniformly distributed over the membrane area (sections with equal membrane area had the same expected number of synapses), and synaptic weights were tuned so that the average single-neuron action potential (AP) firing rate was about five APs per second (this is within the range of typical firing frequencies observed for cortical neurons [63]). As illustrated in Fig 1, a piece of tissue was subdivided vertically into 15 depth intervals (here referred to as ECS subvolumes), which we could picture as spanning from the bottom to the top layer of a layered structure such as cortex or hippocampus. The neurons were positioned so that they occupied the 13 interior subvolumes. The output from all neural segments contained in a specific subvolume were summed, and this gave the total output into the given subvolume. In the neuronal output signal we kept separate track of the different kinds of transmembrane currents, including (i) the net Na+ current, (ii) the net K+ current, (iii) the net Ca2+ current, (iv) non-specific ionic currents, and (v) the capacitive current. For simplicity, we assumed that all unspecified ionic currents in the model [62] (such as leakage currents, synaptic currents, and currents through non-specific active ion channels) were carried by a single, non-specified anion species X-. We chose to use an anion, because many of the non-specified currents are likely to be mediated largely by Cl- (for further comments on this choice, see Methods and Discussion). The output from the neural population into three selected ECS subvolumes is shown in Fig 2 for the first seven seconds of the simulation. For example, Fig 2A shows the currents into the subvolume (n = 3) containing the somata. Here, we clearly see the brief Na+ (Fig 2A1) and K+ (Fig 2A2) current pulses associated with neuronal AP firing. The current amplitudes were about -30 nA (inward, depolarizing current) for Na+ and 30 nA (outward, repolarizing current) for K+. Generally, the subvolume containing the somata received a higher influx/efflux of ions (Fig 2A) compared to the subvolumes containing the apical trunk (Fig 2B) and apical branches (Fig 2C). These differences have two explanations: First, the somata subvolume contained a larger proportion of the total neuronal membrane area, which generally enhanced the ionic exchange in this subvolume. (Similarly, currents are larger in Fig 2C compared to Fig 2B because the subvolume where the apical dendrites branched out contained a larger membrane area than subvolumes containing a part of the apical dendritic trunk.) Secondly, the somata also had a higher density of Na+ and K+ channels than the dendrites. Accordingly, almost all exchange of Na+ and K+ between the neurons and the ECS occurred in the soma subvolume (compare somatic output in Fig 2A1 and 2A2 to the total neuronal output current in Fig 2D1 and 2D2). For the other ions (Ca2+ and X-), the dendrites contributed with a larger proportion of the total output. As we just saw, the neurodynamics fluctuated vividly on the millisecond time scale. However, the input statistics was the same throughout the simulation, so that the slow time-scale neurodynamics was essentially stationary (see Methods). To illustrate this, we split the seven seconds of neural simulations shown in Fig 2 into five 1. 4 second time intervals, and averaged the total transmembrane current (IM) over the five respective intervals. Fig 3 shows how the (temporally averaged) transmembrane sources were distributed across tissue depth. The spatial profile of IM was essentially independent of which 1. 4 second interval of activity it was averaged over. The main current source (positive transmembrane current, i. e. , net positive charge leaving the neurons) was found in the soma subvolume (n = 3). The main current sinks (negative transmembrane current, i. e. , net positive charge entering the neurons) were found in subvolumes containing proximal apical dendrites (n = 5,6) and distal, branching apical dendrites (n = 12,13,14). We note that the transmembrane current profile summed to zero across depth, meaning that the sinks and sources balanced each others out (no neuron can be a net current sink nor source). Of course, the neurodynamics and source/sink configurations seen in Figs 2 and 3 depended in a complex way on the particular neuronal morphology and the subcellular distribution of membrane mechanisms and synapses used in the simulations. The main objective of this work was, however, not to analyze these dependencies, but rather to explore how the ECS potential surrounding the neuronal population depended on whether diffusion was included in the simulations of the ECS dynamics. We investigated this for the particular scenario summarized in Figs 2 and 3, which was used in all simulations shown in the following, but with 84 seconds of simulated neurodynamics, and not only the seven seconds depicted in the figures. We note that the neuron model by Hay et al. exhibited a rich repertoire of firing properties, including the occasional dendritic Ca2+ spikes seen in Fig 2B3. We refer to the original work for further details on the model properties [62]. In the following, the focus will be on how the simulated ECS potential (surrounding this given system) depend on whether ECS diffusion is accounted for. Knowing the neuronal output to each ECS subvolume, we used the KNP-formalism to compute the resulting dynamics of ionic concentrations and the electrical potential in the ECS. Typically, ECS potentials are thought to mainly originate from various transmembrane current sources [16,17]. Here, we explored whether diffusive currents in the ECS could constitute an additional source. Fig 4A1–4A4 illustrates the dynamics of the ECS potential in two selected subvolumes (soma, n = 3, solid line; apical dendrite, n = 13, dashed line) due to the neuronal activity shown in Fig 2. Similarly, Fig 4B and 4C show the field currents and diffusive currents (respectively) from subvolume n = 3 to n = 4 (solid line) and from subvolume n = 13 to n = 14 (dashed line). For simplicity, we in the following discussion refer to the current from n = 3 to n = 4 as the current out from the soma subvolume, and the current from n = 13 to n = 14 as the current out from the apical dendrite subvolume. The first column (1) of Fig 4 shows the time course of these variables over the full simulation, while the remaining columns (2–4) show the time course over selected, shorter (40 ms) time intervals, which include only a few neuronal APs. Red curves represent the scenario without diffusion in the ECS simulations, while blue curves represent the scenario with ECS diffusion included. When we explore the extracellular AP signatures (panels A2–4), we see that they had the same time course as the field currents (panels B2–4), while diffusive currents varied little at this fast time scale (panels C2–4). Diffusive currents thus had no impact on the fast temporal dynamics, and the AP signatures resembled those previously studied in models based on volume-conductor theory, where diffusive currents are neglected [34]. Somatic AP generation was due to an inward (depolarizing) current into the neuron followed by an outward (repolarizing current). Since the sum of transmembrane currents over the neuron as a whole (all ionic + capacitive currents) must be zero at all times, the dendritic branches experienced the opposite current configuration during the APs (outward currents followed by inward currents). Therefore, AP signatures in the apical ECS subvolume (dashed lines in Fig 4A2–4A4) had the opposite temporal profiles compared to what we observed in the soma subvolume (solid lines in panels A2–4). Although the AP signatures were of the same order of magnitude in the soma and apical subvolumes, ECS field currents out of the soma subvolume were generally much larger than field currents between neighboring dendritic subvolumes (panels B2–4). The explanation lies in the spatial distribution of transmembrane inward and outward currents, and the rather unique role played by the soma. For example, a local inward current to the soma returned to the ECS in a widespread manner, i. e. , it was distributed over the entire dendritic tree. Neighboring dendritic subvolumes therefore had similar AP signatures, implying that the ECS voltage differences (and therefore the field currents) between them were small. The diffusive currents varied at a much slower time scale compared to field currents (Fig 4C). This was due to the slow time scale at which ion concentrations varied (as we shall explore further below). The diffusive current out of the soma region reached a peak value after around 30 seconds, after which it decreased slowly. The concentration build-up was slower in the subvolumes containing apical dendrites, and diffusive currents were smaller there, and still increasing at the end of the simulation (panel C1). In the early part of the simulation, when diffusive currents were small, the ECS potential V was close to identical in the cases with and without diffusion (panel A2). However, as diffusive currents built up, they did have an effect on V, which was shifted to more negative values in the simulation with diffusion included compared to case without ECS diffusion (panel A3–A4). Towards the end of the simulation, diffusion had shifted V by about -0. 2 mV in the soma subvolume. In the following, we shall explore this process in further detail. Diffusive currents in the ECS are proportional to concentration gradients in the ECS. To gain insight in the slow dynamics of the diffusive currents, we must therefore investigate the ECS ion-concentration dynamics. In our simulations, ECS concentrations varied due to ionic output from the neurons. Fig 5 shows how the ECS concentration varied over the tissue depth at selected time points. The deviations from the initial concentrations became gradually larger throughout the 84 second simulation, illustrating the slow time scale of ion-concentration dynamics in the ECS. When diffusion was not included in the ECS simulations (Fig 5A), ionic transports were solely due to electrical migration, and were not biased towards following concentration gradients of distinct species. In this case, the ECS concentration profiles predominantly reflected the distribution of neuronal sources. For example, somatic AP generation caused a sharp decrease in the Na+ concentration and a corresponding increase in the K+ concentration in the soma subvolume, while the Na+ and K+ concentration changes were relatively small outside this subvolume (Fig 5A2 and 5A3). We note that the ion-concentration changes in the soma subvolume were unphysiologically high in the no-diffusion case. However, this was of no concern in the current study, since ion concentrations had negligible impact on the ECS dynamics in the case where diffusion was not included. (In this case the only effect on the ECS potentials came from the concentration dependence of the ECS conductivity, see Methods, Eq 11. However, for the present case the conductivity changes were found to be too small to have a visible impact on V in the simulations, see Discussion). With diffusion included in the ECS simulation, the ion-concentration gradients across the depth of the piece of tissue became smoother (Fig 5B). For example, a fraction of the K+ expelled during somatic AP firing diffused out of the soma subvolume, and distributed across the entire tissue volume. In this case, the K+ concentration in the soma subvolume increased from a baseline level of 3 mM to slightly above 10 mM during the 84 second simulation, accompanied by a similar reduction in the Na+ concentration. These concentration shifts were within the range that can be expected under non-pathological physiological conditions (for K+, the limiting concentration between non-pathological and pathological conditions is typically estimated to be between 10 and 12 mM [7]). The buildup of ECS concentration gradients explains the temporal development of the diffusive current that we observed in Fig 4C1. Early in the simulation, the diffusive current out of the soma subvolume (i. e. , from n = 3 to n = 4) increased in an approximately linear fashion with time. This was because the local ion concentration in the soma subvolume (n = 3) increased in an approximately linear fashion due to the high neuronal output/input in this subvolume. As the ion-concentration gradients built up, diffusion from n = 3 to n = 4 increased, and the concentration increase in n = 3 became sublinear. Eventually, diffusion tended to smoothen out the ECS ion-concentration gradients (Fig 5A), and after about 30 s, diffusion between n = 3 and n = 4 experienced a slight decrease. A similar process took place over the entire tissue depth, but was slower further away from the soma, as the transmembrane ionic exchange was smaller there. In the apical dendrites (i. e. , diffusion from n = 13 to n = 14), the diffusive current still increased in a close to linear fashion at the end of the 84 second simulation (Fig 4C1). Due to the slow nature of diffusive currents, we proceeded to investigate the slow time scale dynamics of the ECS potential. To do this, we took the time series of V (plotted for selected subvolumes in Fig 4), and split it up in five equal time intervals of 16. 8 second duration (adding up to the total simulation time of 84 seconds). Next, we took the temporal average of V in these five intervals and obtained a (very) low-pass filtered version of the ECS potential. The results are displayed in Fig 6 showing how the low-pass filtered V was distributed across the tissue depth in the cases without (Fig 6A2) and with (Fig 6B2) diffusion included in the ECS simulations. We first investigate the ECS voltage gradients obtained in the case where ECS diffusion was not included in the simulations (Fig 6B). In this case, there was an ECS voltage drop (of about 1. 3 mV) from the soma subvolume to the subvolumes containing the apical dendrites. The drop in V was consistent with the neuronal source/sinks configurations that we observed earlier (Fig 3): Since the main neuronal current source (transmembrane current entering the ECS) was found in the soma subvolume (n = 3), while the sinks (transmembrane current leaving the ECS) were located higher up along the apical dendrites, there had to be an ECS current in the positive z-direction (corresponding to a negative voltage gradient in this direction) to close the current loop between the sources and sinks. Similar V profiles have been seen experimentally where sustained voltage profiles which vary by a up to several mV at spatial scales of millimeters have been seen in cortex [1,3], hippocampus [26] and in the spinal cord [23]. We also note that the neuronal current sources/sinks were effectively constant at this slow time scale (Fig 3), meaning that they were essentially the same in all the five different time intervals in Fig 6. We would then a priori expect the ECS current to be constant over time as well. Without extracellular diffusion, this would in turn imply that also the ECS voltage gradient should remain constant throughout the simulation, which is indeed what is observed in Fig 6A2 (lines are on top of each others). With diffusion included in the ECS simulations, the situation became more complex (Fig 6C). The gross features of the ECS voltage gradient resembled what we saw in Fig 6B. The similarity was not surprising, since the neuronal sources were identical in the two cases. However, with diffusion included, the ECS potential gradients no longer remained constant throughout the simulation (Fig 6C). The time-dependent variations were most pronounced in the soma subvolume where the ECS potential decreased by about 0. 2 mV over the time course of the simulation. This shift in V was caused by diffusive currents along the ion-concentration gradients that built up during the simulation, and was the same shift that we previously observed in Fig 4A4. A detailed physical interpretation of the diffusion-induced shifts in the ECS potential is provided in the following subsection. To obtain a more thorough understanding of the interplay between the potential V and diffusive currents, we next plotted the ECS fluxes of all ion species (K+, Na+, Ca2+, and X-) in the cases without and with extracellular diffusion (Fig 7). Also here, the focus was on the long time-scale dynamics, and we compared the time-averaged fluxes taken over five 16. 8 second time intervals (same procedure as used for V in Fig 6). In the rightmost column in Fig 7, we have also plotted the total electrical ECS current associated with the ionic fluxes (the definition is given in the caption of Fig 7). When ECS diffusion was not included in the simulations, all ion transport in the ECS were due to the electrical field (Fig 7A). In that case, most of the transports were mediated by the most abundant ion species in the ECS, which in our simulation were Na+ and X-. Due to the negative potential gradient between the subvolumes containing the soma and apical dendrites (Fig 6C), the positively charged Na+ ions were driven away from the soma subvolume, while the negatively charged X- ions were driven towards the soma subvolume. Both these ion fluxes amounted to a net electrical current away from the soma subvolume, i. e. , a positive current in subvolumes above the somata (n > 3) and a negative current in subvolumes below the somata (n < 3). In simulations including extracellular diffusion we plotted the ECS flux densities due the electrical field (jf) and diffusion (jd) separately (Fig 7B1 and 7B2), as well the total flux density (jf + jd, Fig 7B3). As AP firing evoked a decrease/increase of Na+/K+ in the soma subvolume, ECS diffusion drove Na+ into this subvolume, while it drove K+ out of this subvolume (Fig 7B2). As these two cation fluxes were oppositely directed, the net diffusive charge transport (id/F) was smaller than the charge transported by Na+ and K+ separately. However, the diffusive fluxes still gave rise to a net electrical transport of the same order of magnitude as the field-driven current, especially around the soma subvolume (compare current densities in panels B1 and B2 in Fig 7). The ionic fluxes in the ECS differed quite significantly between the cases with and without ECS diffusion (compare flux densities in panels A with B3 in Fig 7). However, the net electrical current in the system were identical in the two cases (compare current densities in panels A and B3). This can be understood from basic electric circuit theory: As the neuronal transmembrane sources/sinks were identical in the two cases, the same had to hold for the net extracellular current. Otherwise, the current loop would not be completed. This leads to the following key insight: Since the net electrical current density (itot = if + id) was independent of whether diffusion was present in the model or not, an increase in id had to be accompanied by a corresponding decrease in if, and vice versa. A time-dependent variation of diffusive currents therefore by necessity evoked a time-dependent variation of the field currents (Fig 7B1) ). As if was proportional to the voltage gradient, this in turn implied that the ECS voltage gradients varied with time, as observed in Fig 6. So far, we have demonstrated that diffusive currents can have quite substantial effects on ECS potentials, at least on a slow time scale. As a next inquiry, we would like to know the frequency range in which diffusion can be expected to have an effect on recorded ECS potentials, and in particular whether diffusion can be expected to affect experimental LFP recordings where the low-frequency cut-off typically ranges from 0. 1 Hz to 1 Hz (see e. g. , [29,30]). We limited this study to ECS potentials recorded in the soma subvolume, where the diffusive effects were most pronounced in our model. Fig 8 shows the power spectral densities (PSDs) of the ECS potential recorded outside the somata (n = 3), where V was obtained as in the above simulations in Figs 4–7). To explore the development of the PSDs over the time course of our simulation, we split the 84 second time series of V into four 21 second intervals, and computed the PSD for these time intervals separately. A first observation is that the PSDs for the simulations without (red lines) and with (blue lines) ECS diffusion differed dramatically for the lowest frequencies, where the presence of diffusion boosted the PSD by up to several orders of magnitude. Contrarily, for the highest frequency components the PSDs were close to identical in the cases without and with diffusion (red and blue lines overlap). This was as expected from our previous analysis where we saw that diffusion was important for the slow, but not the fast system dynamics (Fig 4). The cross-over frequency for which the diffusion contributed negligibly to the PSD, was for all four time intervals depicted in Fig 8 seen to be in the frequency range between 1 and 10 Hz. Extracellular diffusion was thus found to have effect on the PSD for frequency components well within the range typically considered in recordings of LFPs in vivo [29,30]. The PSDs obtained with no ECS diffusion (red lines) were quite constant throughout the simulation, while the PSDs obtained with ECS diffusion included (blue lines) were generally higher for the earliest time intervals (compare panels A and D). To provide a hand-waving explanation to the latter, we start by noting that the contribution of diffusion to the local PSD essentially depended on the absolute value of the temporal variation of local ion concentration (i. e. , on | c ˙ k |, see S1 Appendix), which in turn depended on two competing processes. The first process was the local neuronal output of ion species k, which was roughly constant at the long timescale considered here. The second process was ECS transportation of ion species k out from/into the local region. Generally, these two processes had opposing effects on the local ion-concentration dynamics (i. e. , when neurons expelled K+ into a given subvolume, ECS transports tended to drive K+ out from that subvolume). Early in the simulation, ECS concentration gradients (and thus ECS diffusive transports) were small, and the time development of the local concentration was approximately proportional to the neuronal output. At a later stage, ECS concentration gradients had built up, and the competing diffusive process had increased. Then local concentrations changed more slowly with time. In the rather complex scenario studied so far, transmembrane and extracellular currents interacted (as is, of course, the case in real brain tissue). However, diffusive fluxes and currents in the ECS can in principle exist even without on-going neuronal sources, provided that there are concentration gradients present in the ECS. To improve our understanding of diffusion-generated potentials, we explored them also in such a simplified scenario. For simplicity, we used the same simulation as above (Figs 2–7) to generate reasonable ECS concentration gradients needed in the simplified scenario. However, this time we turned off the neuronal current sources midways in the simulation (i. e. , after 42 seconds), and analyzed the ECS dynamics in last 42 seconds of the simulation when the ECS dynamics was solely due to diffusion along the concentration gradients that had built up during the first 42 seconds of the simulation. For this scenario, only the simulations with ECS diffusion included gave non-trivial results (when extracellular diffusion was not included, the ECS voltage gradient instantly turned to zero when the neuronal current sources were removed, and the extracellular ion fluxes immediately stopped). This can be easily understood from the current conservation laws upon which the KNP formalism was based, stating that the sum of currents into an ECS compartment should be zero (Fig 1C). In the simplified scenario, there were no transmembrane sources after 42 seconds, and with no diffusive currents between ECS subvolumes, the field currents (and thus voltage differences) between ECS subvolumes must by necessity also be zero. The simulations with ECS diffusion included are shown in Fig 9. Panels A2–5 show the ECS concentration profiles at selected time points after the neuronal sources were turned off at t = 42 s. Initially (i. e. , at t = 42 s), ionic concentrations of Na+, K+, Ca2+ and X- in the soma subvolume had been shifted by approximately -5. 1 mM, 6. 0 mM, -0. 1 mM and 0. 7 mM, respectively, relative to the baseline concentrations. We note that these shifts fulfilled the requirement of local electroneutrality, i. e. , did not correspond to any net change in the local charge density: Σk (zkck) = (−5. 1 + 6. 0 − 2 × 0. 1 − 0. 7) mM = 0. Here zk and ck are the valence and concentration, respectively, of ion species k. The deviations from baseline concentrations were smaller outside the soma subvolume, and the concentration gradients out of the soma subvolume were quite steep. Diffusive currents along these gradients gave rise to a diffusion potential, which at t=42 s peaked in the soma subvolume where V was about -0. 17 mV (Fig 9B2). Diffusion-evoked voltage gradients like this are well understood, and have been observed in many systems with spatial variation in ion composition [3,20–22]. With no neuronal sources present, the ECS concentration gradients were gradually smoothed over time (i. e. , for t> 42 s). Consequently, the ECS voltage gradients decayed. At the end of the simulation (i. e. , for t = 84 s), V was about -0. 05 mV in the soma subvolume. The PSD corresponding to this decay process is depicted by the black lines in Fig 9B3 and 9B4. Since the concentration gradients became gradually smoother, the power was generally higher during the first 21 s after the neurons were turned off (Fig 9B3) than in the proceeding 21 s (Fig 9B4). In both cases, the PSDs were very close to a 1/f2 power law (the fitted power-law coefficients were 1. 998 in panel B3 and 2. 02 in B4). This so-called Brownian-noise power law essentially follows from an exponential decay of local ion concentrations, and can be derived analytically (see S1 Appendix). For comparison, we also show the PSDs of the simulation with neuronal sources included (the red and blue lines in Fig 9B3 and 9B4 are the same as in Fig 8C and 8D, respectively). Also in the presence of neuronal sources, the electrodiffusive ECS process roughly followed a 1/f2 power law for low frequencies where diffusion dominated (blue line and black line close to parallel for f < 10 Hz). Comparing the blue and black lines, we further note that the removal of neuronal current sources at t = 42 s increased the low-frequency components of V, especially during the first 21 second time interval after the time of the sources offset (Fig 9B3). To explain this, we may recall that the diffusive power spectrum is proportional to the absolute value of the temporal variation of local ion concentration (| c ˙ k |). As argued above, this value depends on the balance between two competing processes, i. e. , the local neuronal output of ion species k and the ECS transports of ion species k out from/into the local region. The observation in Fig 9B3 simply implies that the local concentration approached the baseline levels faster when the neuronal sources were turned off (black line) than it diverged from the baseline level in the case when the neuronal sources were kept on (blue line). In reality, transmembrane current sources and ECS transport processes do interact, and the correct electrodiffusive PSD is predicted by the blue line in Fig 9B3 and 9B4. Likewise, the predicted maximum frequency that will be affected by diffusion is in the frequency range 1–10 Hz where the red and blue lines in Fig 9B3 and 9B4 merge. However, we still believe that the study of the simplified decay process process (with neuronal sources turned off) provide useful insights to how ECS diffusion can affect the PSD. Firstly, the simplified ‘decoupled’ model nicely illustrated that ECS diffusion gave rise to a 1/f2 contribution to the PSD, as we saw above. Secondly, we propose that the crossing point between the PSD obtained for the diffusive process alone (with concentration gradients representable for what one typically see in the system) and the PSD obtained from neurodynamics when diffusion was not included (black vs. red line in Fig 9B3) may serve as a crude estimate of the maximum frequencies for which diffusion can be expected to influence the PSD. For example, the crossing point between the red and black line in Fig 9B3 was found in the frequency range 1–10 Hz, which agreed with the frequency range where the blue and red lines merged. We will provide further arguments for the usefulness of the simplified scenario further in the Discussion. So far, we conclude that in the current model, diffusive processes affected ECS potentials for frequencies up to several hertz.
The simplified model set-up used here have several limitations. Firstly, V was computed as an averaged value over a large ECS volume, and comparison between this and experimental recordings of V with point-like electrodes with small contacts is not straightforward. Secondly, brain tissue contains many types of neurons, which are distributed with somata in different depth layers (see, e. g. , [18,69,78]), whereas we only included one. Thirdly, we did not include synaptic connections between neurons. Such connections could induce a level of synchrony in the neuronal firing, which likely would influence the power spectrum of the ECS potential [16,35]. Fourthly, we assumed that spatial variations in the electric potential and ion concentrations occurred only in one spatial dimension. This is clearly not strictly true, and some aspects of the estimated power spectra are likely to depend on the three-dimensional nature of the real system. Fifthly, the presently used multicompartmental neuronal model [62] (together with most other available multicompartmental models) does not include ionic uptake mechanisms such as Na+/K+-pumps. Such mechanisms, along with glial uptake mechanisms [28,79], would generally act to maintain the ECS ion concentrations closer to the baseline levels than what we predicted with our model. These shortcomings are discussed in further detail below. In the current section we discuss the predictions that we made regarding diffusion-generated electric potentials, and to which degree these can be expected to reflect realistic experimental scenarios. To clarify the discussion, we start by labeling the three situations that we have studied P1, P2, and P3, respectively. We refer to Fig 9B3 and 9B4, where all three situations (P1–P3) are represented. The blue line (P1) represents ECS dynamics surrounding an active neuronal population in the realistic scenario described by the full electrodiffusive formalism. The red line (P2) represents the situation where ECS diffusion was neglected so that the ECS potential was given exclusively by the distribution of transmembrane sources (cf. , standard volume-conductor theory [32]). Finally, the black line (P3) represents the situation where the neuronal sources had been turned off, so that the ECS potential was driven exclusively by concentration gradients in the ECS. The concentration gradients could in principle be imposed as an initial condition in the system, independent of the neural model, meaning that P3 and P2 were essentially independent processes. As we shall see below, this independence is useful for analyzing our results, and for comparing them to previous studies. A common starting point for the estimation of the current-source density CSD (x, y, x) from the ECS potential V (x, y, z) is [37,39]: ∇ (σ ∇ V) = - C S D (1) The left hand side is the divergence of the ECS currents, and an implicit assumption in this equation is that only electric currents driven by the electrical field is present in the ECS, i. e. , solely Ohmic current densities given by if = −σ∇V. If also diffusive ECS currents were accounted for, the corresponding equation would be: ∇ (σ ∇ V) - ∇ i d = - C S D, (2) where the diffusive current density is a function of ionic concentrations in the ECS (see Methods, Eq 9). The use of Eq 1 for predicting the CSD could thus lead to a misinterpretation of diffusive ECS currents (if present in the real system) as neuronal current sources. An example where experimental recordings seems to disagree with the standard CSD theory (Eq 1) was reported recently by Riera et al. [30] who found that the estimated instantaneous current-source density (CSD) from recorded ECS potentials did not sum to zero over the volume of the barrel column. According to the standard CSD-theory, this would indicate the presence of a non-zero current-source monopole on a mesoscopic (cell population) scale. The possible origin of these apparent current monopoles was later debated [30,103–106]. A non-negligible diffusive source term, cf. Eq 2 could be one (of several) possible explanations of this discrepancy between experiments and original CSD-theory. The model presented here was a simplified one, both in terms of using a 1D geometry and in terms of neglecting several neuronal and glial mechanisms that would likely contribute to the generation of the LFP. A future ambition is to expand this framework to a 3D model that also accounts for more of the complexity of neuronal tissue, and includes effects of neuronal and glial ionic uptake mechanisms (ion pumps). A 3D version of the KNP framework could ideally be combined with existing, comprehensive simulators of large neuronal networks such as the Blue Brain simulator [69]. We believe that such a framework would be very important for the field of neuroscience as it not only would be useful for exploring how diffusive currents can have an impact on ECS potentials, but also to simulate various pathological conditions related to ion-concentration dynamics in neural tissue [9,13–15,101,107].
What we have here coined the Kirchhoff-Nernst-Planck formalism, was originally developed for computing the intra- and extracellular dynamics of ion concentrations and the electrical potential during astrocytic K+ buffering [28]. In the current application, it was only applied in the ECS (the intracellular space was handled with the NEURON-simulator). For simplicity, we assumed that spatial variation only occurred in one spatial direction (z-direction), and thus that we had lateral homogeneity of all state variables. In the current work, the KNP formalism was used to predict the extracellular ion-concentration dynamics and electrical potential surrounding a small population of ten pyramidal cells. The neural simulation used in this study was briefly introduced in the Results section, but is presented in further detail here. Simulations on the pyramidal cell model by [62] was run the NEURON/Python simulation environment [60]. The ECS dynamics was computed separately with the KNP formalism, using the neuronal output/input as an external input time series. The KNP model was implemented in MATLAB (http: //se. mathworks. com/). The MATLAB code (along with the neuronal input time series) will be made publicly available at ModelDB (http: //senselab. med. yale. edu/modeldb).
When electrical potentials are measured in the extracellular space (ECS) of the brain, they are interpreted as a signature of neural signalling. The relationship between the ECS potentials and the underlying neuronal processes is often studied with the aid of computer models. The ECS potential is typically assumed not to be affected by diffusive currents in the ECS, and existing models therefore neglect diffusion. However, there may be scenarios where this assumption does not hold. Here, we present a new computational model which explicitly models ion-concentration dynamics in the ECS surrounding a neural population, and which allows us to quantify the effect that diffusive currents have on the ECS potential. Using this model, we simulate a scenario where a population of pyramidal neurons is active over a long time, and produces large, but realistic concentration gradients in the ECS. In this scenario, diffusive currents are found to influence the ECS potential at frequency components as high as ten hertz. Unlike previously believed, we thus predict that there are scenarios where recorded local field potentials (LFPs) are likely to contain signatures not only of neural activity, but also of ECS diffusion.
medicine and health sciences
nervous system
membrane potential
ions
electrophysiology
neuroscience
simulation and modeling
computational neuroscience
neuronal dendrites
physical chemistry
research and analysis methods
animal cells
chemistry
cellular neuroscience
neuronal morphology
cell biology
anatomy
synapses
physiology
neurons
single neuron function
biology and life sciences
cellular types
physical sciences
computational biology
neurophysiology
Effect of Ionic Diffusion on Extracellular Potentials in Neural Tissue