Relationship to other modelling schemes.

What we here have coined the KNP formalism was developed in previous work where we derived a mathematical formalism for simulating buffering of extracellular K⁺ by astrocytes [28, 57]. A very similar formalism was developed (in parallel) in the heart-cell community in the context of a model of ischemia [58]. Our buffering model accounted for electrodiffusive processes in the intra- and extracellular domain, and was essentially an expansion of the previous model by Qian and Sejnowski [64], which only considered the intracellular domain.

The KNP formalism represents a simplification of the computationally expensive PNP-solvers, which derive the local potential from Poisson’s equation (e.g., [47–51]). The PNP system has been thoroughly analyzed in a series of previous works by Mori, who also proposed a series of simplified models and studied their validity under different conditions [53–56, 65]. The KNP formalism can be regarded as a simplified version of the electroneutral model proposed by Mori [52, 54, 65]. In the current application of the KNP model, we tailored the formalism to study transport processes in neuronal tissue at a relatively large spatiotemporal scale.

In the original application of the KNP formalism, we modelled both the intra- and extracellular space explicitly [28]. In the current application, we introduced a hybrid modelling framework, where the KNP formalism was only applied to the ECS, while the intracellular dynamics was computed with the NEURON simulator [59, 60]. The NEURON simulator is an efficient standard tool for computing the dynamics of morphologically detailed neurons. It can be combined with algorithms for handling the intracellular dynamics of selected ion species due to local transmembrane influxes/effluxes and decay processes. (In most neural models, such algorithms are typically an exception used only for the signalling molecule Ca²⁺, see e.g., [62, 66–68].) However, intracellular electrodiffusive processes are so far not an integral part of the NEURON simulator. A limitation with the hybrid scheme is therefore that electrodiffusive processes are only accounted for in the ECS domain, where the KNP-formalism is used. An important advantage with the hybrid scheme is that it lends itself to be used as a supplement to compute the ECS dynamics (of ion concentrations and the electrical potential) of the multitude of already available neural or neural network models based on the NEURON simulator (such as, e.g., the Blue Brain simulator [69]). Previous tools developed to compute extracellular potentials from NEURON based models such as LFPy [61] have not incorporated effects of ionic diffusion in ECS.

Electroneutrality assumption.

Along with other electroneutral models [53–56, 58, 65], the KNP formalism provides a means of deriving the local potential V from the physical constraint that the (intra- and extracellular) bulk solution is electroneutral [28]. This approximation was used as early as in 1890 by Planck, who described electrodiffusion in electrolytes [70]. In the current application, this approximation means that any nonzero local charge density within the system is identified as a charge that sits on a capacitive membrane and uniquely determines the local transmembrane potential of an excitable cell. Put differently, the KNP scheme assures that the sum of ionic currents into a given tissue sub-volume equals the sum of capacitive (non-ionic) currents over the cellular membranes that populate the sub-volume (as illustrated in Fig 1C), so that no net charge is found in the bulk.

The assumption that bulk solutions is electroneutral is not strictly true, as has been the topic of many discussions (see e.g., [45, 71]). Indeed, Fig 9 showed that in the presence of diffusion, we could obtain a nonzero voltage gradient even in the absence of neuronal sources, an observation which is incompatible with the notion of a strictly electroneutral ECS. However, it has been shown that invoking the electroneutrality assumption is equivalent to invoking the limit of the exact PNP treatment when charge-density-dependent effects become small [71], and that the electroneutral model works as an excellent approximation at spatiotemporal scales larger than microseconds and micrometers [52, 65].

Diffusion potentials.

Diffusion-generated potentials are well known in electrolyte theory. Often they are referred to as liquid junction potentials, since they are most pronounced at the boundary between two solutions of different ion composition [20–22, 72].

In reality, the genesis of liquid junction potentials is a three-step process that requires (i) initial ion-concentration gradients that are such that diffusion will drive a net electrical charge in some direction, (ii) a charge separation associated with the diffusive process, and (iii) an electrical potential that arises from the charge-separation process, and opposes further charge separation. This diffusion-generated potential (iii) represents a quasi steady-state scenario where electrical drift and diffusive drift are opposing and in equilibrium. Simplified equations for computing diffusion potentials include the Henderson equation and the Goldman-Hodgkin-Katz equation (see e.g., [22], and S2 Appendix). The relaxation towards this quasi-steady state occurs very rapidly, i.e., on the nanosecond timescale [73]. Furthermore, the number of ions that constitute the net charge density during equilibrium is about nine orders of magnitudes smaller than the number of ions present [72]. The KNP formalism bypasses the rapid equilibration process by assuming that the quasi-steady state is reached instantaneously, and derives the value for V associated with the equilibrium state. In doing so, the KNP formalism implicitly neglects the tiny local charge separation associated with the charge relaxation process.

To get an intuitive understanding of the diffusion potential, it may help to compare it with the (in neuroscience) more familiar cellular resting potential, which is typically computed from the Goldman-Hodgkin-Katz (GHK) equation (see S2 Appendix). The GHK-equation predicts the equilibrium potential between two compartments A and B with different ion compositions, i.e., the potential difference at which the diffusive current from A to B equals the field driven current from B to A. When the resting potential is computed, the compartments A and B represent the inside and outside of an excitable cell membrane. However, if we let A denote the subvolume n = 1 where the ion concentrations (by constraint) had the baseline values (Table 1), and let B denote the compartment n = 3 with concentrations as in Fig 9 at t = 42 s, the GHK equation (see S2 Appendix) predicts a potential difference of −0.17 mV between the two compartments. This agrees with the value we got in Fig 9 when neuronal sources were turned off at t = 42 s.

Extracellular conductivity.

Unlike in volume-conductor theory, where a fixed value typically is used for the ECS conductivity [18, 32–34], the KNP formalism models the conductivity (σ) as a function of the number of free ionic charge carriers, weighted by their mobility and valence (cf., Eq 11). In neural tissue, the main charge carriers are typically believed to be K⁺, Na⁺, and Cl^-. The model included K⁺, Na⁺, Ca²⁺ and an unspecified anion species (X^-). The latter essentially represented Cl^- in the biological system and was given the same baseline ECS concentration and diffusion constant as Cl^-. The model thus included the main charge carriers, and with the initial ion concentrations that we used, we obtained an ECS conductivity (σ = 0.76 S/m).

In the literature, there are quite some variations in values that are given for the ECS conductivity, and also variations in how this quantity is defined. In the current study, we used the porous medium approximation [74], and explicitly accounted for the fact that ECS currents only go through a volume fraction of about 0.2 of the tissue volume (see methods). However, it is common to rather define an apparent tissue conductivity, σ′, which is defined using the tissue as a whole as reference volume for ECS currents [4]. In our case, the apparent conductivity was thus σ′ = ασ ≃ 0.15 S/m. For comparison, Chen and Nicholson found an apparent conductivity of σ′ = 0.1 S/m [4], while other computational studies of local field potentials and current-source densities have used values σ′ ∼ 0.3 S/m [18, 32–35, 75]. Our estimate thus lies between the previously estimated values for σ′, and is relatively close to the value used by Chen and Nicholson [4].

Relative variations in ion concentration were quite small in the simulations studied here (at least when it comes to the most abundant species). In addition, such variations tended to be asymmetric (e.g., decreases in K⁺ were accompanied by increases in Na⁺), meaning that variations in the net number of free charge carriers were even smaller than variations in individual ion species. Therefore, σ only varied by a few percent relative to the initial value during simulations. As verified in additional test simulations where σ was pegged at the initial values, these variations had no significant effect on the simulation results.

Albeit its concentration-dependent magnitude, the conductivity (as defined here) was essentially a pure, resistive conductor, i.e., it was independent of the frequency of currents passing through it. A frequency independent conductivity finds support in recent experiments [76] (although there are also experiments that have indicated otherwise [77]).

No feedback from extracellular space to neurons.

To compute the ion-concentration dynamics in the ECS (Fig 5), we counted the number of ions exchanged between the neurons and the ECS in simulations of the multicompartmental neural model [62]. For simplicity, we assumed that that there was no feedback from the ECS dynamics to the neurons. That is, we did not account for changes in neural reversal potentials due to changes in ECS ion concentrations [9, 12], or ephaptic effects of ECS potentials on neuronal membrane potentials [53, 80–82]. Such feedback mechanisms would likely influence the neurodynamics. However, we do not believe this to be a main concern with the current study, as the main focus was not the neurodynamics as such, but rather the ECS dynamics resulting from it. Having no feedback also gave us the advantage that we could have exactly the same neurodynamics when comparing the ECS dynamics in the cases without or with diffusion in the ECS.

Unspecific ion species.

Only a subset of the transmembrane currents in the multicompartmental neuron model [62] were ion specific, and we therefore assumed that all non-specific currents were mediated by an unspecified anion species (X^-). For simplicity, we assumed that all non-specific currents (including the leakage current, synaptic currents, and the currents through non-specific, active ion channels were mediated by the same ion species X^-.

Although this is an inaccurate assumption (e.g., leakage currents are composed of several ion species, and not only one, and the ion channel I_h in the neural model [62] is in reality a cation-channel), it was not critical for the simulation outcome. One reason for this is that the ECS ion concentrations did not have any significant effects on the ECS conductivity (as clarified above), so that the field currents depended little on the composition of ions in the ECS. As for the diffusive currents, the concentration gradients were most dramatic for K⁺ and Na⁺, and the remaining ion species (X^- and Ca²⁺) gave only minor contributions to ECS diffusion. A subdivision of X^- into different ionic species would therefore expectedly not change our qualitative findings.

Simplified neuronal model system.

Albeit internally consistent, the model system represented a crude simplification of the complexity of real tissue, where an intricate circuitry of many different neuron species are likely to contribute to the ECS dynamics. The small population of 10 pyramidal cells used in the current study will likely create a bias towards strong concentration gradients surrounding the soma subvolume (n = 3 in Fig 1). In addition, brain tissue contains neuronal [83] and glial [2, 3, 6] uptake mechanisms (in particular Na⁺/K⁺-exchangers), which were not included in the current model. Such uptake mechanisms work to keep ECS concentrations close to baseline levels, and significant changes in ECS ion concentrations are therefore likely to occur only in cases when the neuronal activity level is too intense for such clearance mechanisms to keep up. The ionic concentration gradients predicted in Fig 5 are therefore likely to be an overestimation of the ion-concentration gradients that would realistically build up during the relatively moderate AP firing activity of the small neuronal population considered here. For comments on how non-included cellular mechanisms could influence the main findings in this work, we point to the discussion below on effects of ECS diffusion currents on measured PSDs.

Volume-averaged potential.

Recorded ECS potentials depend on the distances between the recording electrode and the neuronal current sources [31]. For example, ECS signatures of APs are only large in the vicinity of the neural membrane, while slower signals can have a longer spatial reach [33, 35]. A direct comparison between V as determined by the formalism presented here (averaged over a ECS subvolume), and V measured by point electrodes, would require a generalization of the model to three spatial dimension, using a relatively fine spatial resolution. However, also in the present implementation, V was determined from current conservation laws, and followed the same time course as (Fig 4) the ECS signals seen in previous studies [33, 35]. In addition, the estimates of V in the current work showed sustained ECS profiles (Fig 5) that were qualitatively similar to those observed experimentally [1, 3]. We thus believe that the large scale (volume averaged) V considered in the current study represents a useful quantity for assessing relative contributions of field currents and diffusive currents at a tissue level.

Correlation effects introduced with population size and geometry.

Two important model assumptions could influence the scaling of the PSD observed in this study: the grouping of N neurons into a joint population output I^M into each ECS subvolume, and the assumption of a 1D system geometry.

Regarding the population size, the N neurons in the current model received uncorrelated synaptic input, but with the same (time-averaged) input statistics. This means that the N neurons produced output where the fast components were uncorrelated (e.g., APs were unique for individual neurons) and the slow components were correlated (e.g., the time averaged output was the same for all neurons). Thus, we would expect the high-frequency part of I^M to sum as uncorrelated noise (the amplitude scales roughly like ), and the low-frequency part to sum as correlated signals (the amplitude scales roughly like N). That is, the PSD would be dependent on the system size, so that an up-scaling of the system (increasing population size N and ECS volume by the same factor) would penalize the high-frequency part of the PSD relative to the low-frequency part. Such population effects have been demonstrated in a previous study, where the high frequency part of the ECS potential was found to scale sublinearly with the number of APs elicited in a volume [33].

Similarly, correlation-related effects could also be introduced with the 1D-assumption. In a hypothetical 3D model, the population in Fig 1 would be surrounded by similar populations. In the 1D model, the ECS currents in the lateral directions were by construction zero, which would be equivalent to having zero gradients of the voltage and concentration in the lateral directions. In a (hypothetical) 3D model, this would only occur in the case when neighboring populations were perfectly correlated, which is then an implicit assumption in the 1D model. This could be a good assumption for slow frequency components, but not for the fast components (i.e., the neighboring populations could share input statistics, but not exact AP-spike times). Hence, a transition to a 3D model would likely also penalize the high-frequency part of the PSD relative to the low-frequency part.

Diffusive currents were mainly found to influence the slow components of the PSD, and it is unclear whether the key findings regarding these would be influenced by choice of population size and model dimensionality. This would, however, be a natural topic for future investigations (see below).

Magnitude of diffusion-generated potential (P3).

Figs 6 and 9 showed that the diffusion-generated potential shifts in the ECS developed on a slow time scale. Slowly varying ECS potentials have also been reported in several experimental studies, and may be of the order of several millivolts [1, 3, 13, 23–26]. Generally, the main source of these slow potentials is not believed to be ECS diffusion, but rather glial buffering currents triggered by increases in ECS K⁺ concentrations [3], i.e., on transmembrane current sources that were not included in the computational model studied here. However, based on recorded concentrations differences between different cortical layers during neuronal hyperactivity, Dietzel et al. [3] estimated (using the Henderson equation [84]) that ECS diffusion could contribute to such shifts by maximally 0.4 mV, a finding that they also verified experimentally in a simplified setup. These estimates depended solely on differences in ion-concentration compositions between different cortical regions, and could thus be compared to our scenario P3. For the concentration gradients built up at the time when the neuronal currents were turned off (t = 42 s), we predicted a diffusion generated potential of about 0.17 mV (Fig 9). This was smaller than, but of the same order of magnitude as the maximal shifts estimated in [3].

The magnitude of diffusion potentials depends on spatial variations of ion concentrations, and may be large in non-biological systems (see e.g., [85]). In brain tissue, however, concentration differences are likely more moderate, and diffusion potentials larger than a few tens of a millivolt will probably be rare.

Power spectrum of diffusion-generated potential (P3).

While the magnitude of slow diffusion potentials finds support in previous experimental studies, no previous study has to our knowledge systematically investigated the PSD associated with their temporal development. In the current study we found that diffusion could have an effect on the PSD for frequencies as high as ∼10 Hz. In the discussion following Fig 9B3, we suggested that we could approximate the frequency range where diffusion had an effect with the frequency range where the PDS of the diffusive process alone (P3) had a magnitude that was similar to, or higher than, the PSD predicted from the transmembrane sources alone using volume-conductor theory (P2). We note that this is only an approximation, i.e., without neuronal sources as in process P3, the ion-concentration dynamics and thus ECS diffusion were not identical to that in the full model (P1). The approximation was useful as it allowed us to compare two independent processes (P2 and P3). Whether the predicted frequency range was realistic, can then be boiled down to a question regarding the realism of the PSDs obtained separately for the neuronal model (P2) and with pure ECS diffusion (P3).

The analysis of P3 led to the clear observation that the (undisturbed) diffusion potential followed a 1/f² power law, as we also predicted analytically (see S1 Appendix). The realism of the diffusion-generated ECS potentials (P3) depend predominantly on whether the ion-concentration gradients as such were physiologically realistic, i.e., independently of which underlying neuronal process gave rise to them. Regarding ECS concentrations, most experimental data are available for K⁺, and the simulated K⁺ concentrations were within the range reported experimentally in different systems [1, 3, 4, 6, 8, 24, 86, 87]. For example, the simulated K⁺ gradients in Fig 5B bear some resemblance to K⁺ concentration gradients seen in the experiments by Cordingley and Somjen, where the K⁺ concentration varied with about 4 mM over the depth of cortex (see Fig. 5 in [1]). In line with previous estimates, the simulated shifts in the concentrations of other included ion species were smaller or of the same order of magnitude as for K⁺ [3, 88], and should also be physiologically realistic.

As for the temporal aspect, the simulated ion-concentration variations occurred at a timescale of tens of seconds, which is also similar to what has been seen in experiments (see e.g., [3]), although faster shifts can be induced under specific stimulus conditions [1]. We thus believe that the PSD obtained with the diffusive currents simulated here should be a realistic prediction of what one could observe experimentally under physiological conditions with large ECS concentration gradients.

In the current work, diffusive potentials arose due to large scale concentration gradients in the ECS bulk solutions. We note that these diffusive effects relate to a different phenomenon than the diffusion-evoked 1/f-filtering effects proposed by Bedard and Destexhe (with co-workers). These authors developed a mean-field description of neuronal tissue (comprising both membranes and ECS), and incorporated diffusion-evoked frequency-filtering effects in terms of a complex Warburg impedance [41–43, 89–91]. They argued that the Warburg effects arises due to highly localized diffusion processes in the membrane-near Debye-layers when charge is transferred from the intracellular to extracellular space [90, 91]. Traditionally, the Warburg impedance has been derived for complex interactions at interfaces between electrodes and electrolytes, where the chemical reactions necessary for transferring charge between the electrode surface and electrolyte requires a continuous, diffusion-dependent, reshuffling of local ions [92, 93]. The physical argument to why similar effects should take place close to neuronal membranes is presently unclear, and a complex conductivity, accounting for such possible effects of membrane-near filtering, was not included in the present model. In any case, such putative Warburg-type effects arising close to the membrane and the diffusion potentials evoked by large scale ECS concentration gradients describe different and complementary effects that ionic diffusion could have on the LFP, and are not a priori in contradiction.

Power spectrum of membrane source-generated potential (P2).

The realism of the membrane-current generated ECS potentials (P2), on the other hand, depend on whether neural population model was sufficiently detailed to generate ECS potentials that one would expect under realistic experimental conditions.

It is generally not trivial to construct models that reproduce experimentally recorded PSDs [17]. Previous studies have shown these to be sensitive to the distribution and balance between excitatory and inhibitory synapses on the neuronal membranes [33, 36]. Furthermore, the LFP is also likely to receive contributions from a large fraction of neurons that are not firing APs, but receiving synaptic input, so that they still participate in generating the low powers of the PSD [33]. In addition, different aspects of the LFP has been found to depend on neuronal morphology, subcellular distributions of membrane mechanisms, and the level of synchrony between neighboring neurons [18, 19, 33, 35, 94–96].

In the case of large deviances from baseline ECS concentrations, also glial buffering mechanisms [3, 16] and neuronal uptake mechanisms such as Na⁺/K⁺-exchangers [83] could constitute additional slow membrane currents that could influence the low frequencies of the LFP. Generally, such mechanisms also act to reduce concentration gradients in the ECS. Whereas the concentration gradients seen in the current simulations (not including uptake mechanisms) were realistic, such gradients would probably under most conditions require a higher neural activity level (P2) than in the current model, since the neuronal output would need to out-compete uptake mechanisms in order to generate ECS gradients.

Inhibitory synapses, inactive neurons, ion pumps and glial buffering mechanisms were not included in the model considered here. In a realistic scenario, it is likely that these mechanisms could enhance the membrane-current induced (P2) contribution to the lowest frequencies in the PSD compared to what we predicted in the current model (see below for further discussion on this).

Power spectrum of the full electrodiffusive model (P1).

The PSDs seen in ECS recordings is a highly complex topic. Generally, the frequency scaling is multifactorial, dependent on state (sleep/wake) [97], neural correlation/decorrelation [35, 95], size of active populations [35], morphology of nearby, active neurons [34, 98], specific activity (spiking frequency/synaptic input) of nearby neurons [33], and possibly by frequency filtering within the extracellular medium [41, 97] and diffusion along extracellular concentration gradients, as seen here.

If different causes for frequency scaling are linearly dependent (mathematically, this means that two processes can be expressed as a convolution), it leads to addition of the individual powers (exponents). For example, if an incoming spike train triggers synaptic currents, and both the underlying processes have a 1/f² frequency scaling, the net scaling will be 1/f⁴. As indicated in Fig 9, the large-scale diffusion process considered in the current work is quite independent of other processes, and the power (exponent) from diffusion does not add to the powers (exponents) of the transmembrane sources. In this regard, the large-scale diffusion effects differ from the Warburg-type filtering effects hypothesized in some other works to take place in thin sheaths surrounding neuronal membranes [90, 91]. This implies that large-scale diffusion will will exhibit its characteristic 1/f² frequency scaling, and will be visible in the PSD only in the frequency range where diffusion is the dominant process (when such a range exists). It also means that the presence of diffusion is not in conflict with observing an undisturbed power law generated by other processes (this would only imply that the other process dominates).

As we have argued above, the PSD predicted for the diffusive process (P3) should be physiologically realistic (for large concentration gradients), while the low-frequency components of the membrane-current induced PSD in our model (P2) may be an underestimation of what would be expected in a real system. If additional (slow) membrane mechanisms were included, we would expect the range of frequencies where diffusion dominated the PSD (and the crossing point between the red line (P2) and black line (P3) in Fig 9B3) to be shifted towards lower frequencies (i.e. lower that the 1–10 Hz found here, and in many cases, possibly to frequencies below the cut-off frequency used in LFP recordings). However, given the steepness (1/f²) of the diffusion-generated PSD, a cross-over frequency below which diffusion dominates the PSD, is still likely to occur, especially under conditions where we can expect large extracellular concentration gradients. By prediction, we would then expect to observe a diffusion evoked 1/f² scaling for some low-frequency range of the PSD.

A 1/f² scaling for low frequencies has indeed been observed in cerebral areas of human patients with epilepsy [99, 100], a pathological condition which is strongly associated with dramatic changes in ECS concentrations [101]. A similar scaling was found in slices from rat hippocampus when epileptic-seizure-like events were induced [102]. It should be noted, however, that also processes other than diffusion can give rise to a 1/f² scaling of the PSD [41, 99], and that the 1/f² scaling in one of the cited studies was originally explained by slow-wave state transitions between up/down states [99].

Continuity equation.

The KNP formalism represents a way of solving the continuity equation for the ionic concentrations ( (mol/m³)), and is here derived for the system sketched in Fig 1B. The ECS is subdivided into a number of N = 15 subvolumes of length l_c and cross section area A_c. Using data from cortex (as in [35]), we assume that the average surface area per neuron is about 300 μm², so that the ten neurons used in our simulations occupy a surface area A_c = 3000 μm². The vertical length of a subcompartments is set to l_c = 100 μm, so that the neuron (having a vertical extension of slightly below 1300 μm) occupy the interior 13 subvolumes.

In each subvolume n, the concentrations of all present ion species k are assumed to be known. Ions may enter the subvolume either via (i) transmembrane fluxes from neurons that exchange ions with the subvolume (J^kM), (ii) diffusive fluxes between neighboring subvolumes (J^kd), or (iii) field fluxes between neighboring subvolumes (J^kf). The formalism computes the ECS fluxes (in mol/s), and can be applied together with any selection of neuronal sources. For now, we assume that the transmembrane fluxes J^kM for all ion species as well as the transmembrane capacitive current (which will be relevant below) are known (e.g., determined from a separate simulation using, e.g., the NEURON simulator [60], as we shall return to later). The continuity equation is (in discretized form): (3) where we have used the notation that J_n−1,n denotes the flux from subvolume n − 1 to subvolume n. The extracellular fluxes are described by the Nernst-Planck equations [108]: (4) and (5) where the factor ψ = RT/F is defined in terms of the gas constant (R = 8.314 J/(mol K)), the absolute temperature (T), and Faraday’s constant (F = 96,485C/mol).

We have used the porous medium approximation, characterized by the parameters α and λ [74]. The parameter α represents the fraction of the tissue volume being ECS, and we used the value α = 0.2 [4, 74]. The prefactor αA_cl_c in eq 3 then equals the ECS volume of a subvolume n. The extracellular tortuosity λ represents miscellaneous hindrances to motion through neuronal tissue [4, 74, 109], and gives rise to a reduced effective diffusion constant where D^k is the diffusion constant for ion species k in dilute solvents. We used the value λ = 1.6 [4], and standard values for the diffusion constants [110]: , , = 0.71 × 10⁻⁹m²/s and (Here, X is an unspecified ion species (see below), for which we used the diffusion constant for Cl^-).

We assume that the edge subvolumes (n = 1 and n = N) represent a background where ion concentrations remain constant. The continuity equation then governs the ion-concentration dynamics in all the N − 2 interior subvolumes. If we include a number K of different ion species, the continuity equation (Eq 3) for n = 2, 3, …, N − 1 and k = 1, 2, …, K gives us K(N − 2) conditions for the K(N − 2) ion concentrations in the N − 2 subvolumes where ion concentrations are dynamically changing. However, the continuity equation also includes N state variables for the potential V_n in all subvolumes (including the edges). We thus need N additional constraints to fully specify the system.

Derivation of extracellular potential.

In the following, we derive expressions for the ECS potential (V_n) based on the principle of Kirchhoff’s current law, and the assumption that the bulk solution is electroneutral [28]. To do this, we multiply the continuity equation (Eq 3) by Fz^k, take the sum over all ion species k, and obtain the continuity equation for electrical charge: (6) Here, we have transformed fluxes/concentrations into electrical currents/charge densities by use of the general relations [108]: (7) (8) (9) and (10) where z^k is the valence of ion species k and F is Faraday’s constant. In Eq 10, we also defined the conductivity (units (Ωm)⁻¹)for currents between two subvolumes n − 1 and n as: (11)

At time scales larger than nanoseconds, bulk solutions can be assumed to be electroneutral [111]. In our scheme, bulk electroneutrality implies that any net ionic charge entering an ECS subvolume must be identical to the charge that enters a capacitive neural membrane within this subvolume. This is also an implicit assumption in the cable equation (see, e.g., [28, 64, 108, 112]) upon which the NEURON simulator is based. With this assumption at hand, the continuity equation for charge (Eq 6) becomes useful for us, as it is governed by a constraint that we did not have at the level of ion concentrations (Eq 3). Electroneutrality in the bulk solution implies that the net charge entering an ECS subvolume (the time derivative of q_n in Eq 6) must be identical to the charge which accumulates at the neuronal membrane and gives rise to the neurodynamics. This means that the time derivative of q_n must be equal to the capacitive current that we know from the NEURON simulator: (12) Thus, q_n (in Eq 6) is not an independent state variable, but an entity given from the NEURON simulation (i.e., an input condition to the ECS). With this at hand, we can rewrite Eq 6) on the form: (13) We now see that Eq 13 is simply Kirchhoff’s current law, and states that the net current into an ECS volume n is zero, cf. Fig 1C. If we insert Eq 10 for I^f, Eq 13 becomes: (14)

We note that was defined as the net ionic transmembrane current (Eq 7), and that it does not include the capacitive current. We further note that Eq 14 for a subvolume (n) depends on the voltage levels in the two neighbouring subvolumes (n − 1 and n + 1), and thus only gives us N − 2 conditions, i.e., one for the N − 2 inferior volumes. We need two additional criteria for the edge subvolumes (n = 1 and n = N). As we may chose an arbitrary reference point for the voltage, we may take the first criterion to be: (15)

As the second criterion, we impose a boundary condition stating that no net electrical current is allowed to pass between the subvolumes n = N − 1 and n = N (i.e. no net electrical current enters/leaves the system from/to the constant background). Since there may be a diffusive current between these two subvolumes ( is not constant), this criterion implies that we must define V_N so that the field current is opposite from the diffusive current (). If we insert for I^f (cf., Eq 10), this condition becomes: (16)

The conductivity (σ) and the diffusive currents (I^d) are defined by ionic concentrations in the ECS, whereas we assumed that the neuronal output (I^cap and I^M) was known. Eqs 14–16 thus give us N equations for the N voltage variables V_n. In matrix form, we can write the system of equations (Eqs 14–16) as: (17) where V is a vector containing the potential V_n in all N subvolumes, and b is a vector with N elements given by: (18) The N × N matrix A: (19) is a tridiagonal matrix. The diagonal above the main diagonal is given by: (20) The diagonal below the main diagonal is given by: (21) The main diagonal is given by: (22)

For each time step in the simulation, we can determine V_n by solving the algebraic equation set: (23) where A⁻¹ is the inverse of the matrix A.

When we ran simulations where diffusion was not included, J^d was simply set to zero in the continuity equation (Eq 3), and I^d was set to zero in the equation where the ECS potential is derived (Eq 18).

Initial conditions.

As initial conditions, we assumed that all ECS volumes were at potential V_n = 0. The initial ion concentrations were also identical in all ECS subvolumes. We used c^K0 = 3 mM, c^Na0 = 150 mM, c^Ca0 = 1.4 mM. These ion concentrations are quite typical for cerebrospinal fluid [113]. To obtain an initial charge density of zero in the bulk solution, we computed that the initial concentration for the unspecified anion should be c^X0 = 155.8 mM: With this value, we get that local charge density ρ/F = ∑z^k c^k0 = (1 × 3 + 1 × 150 + 2 × 1.4 − 1 × 155.8) mM = 0. This value for c^X0 is close to typical ECS concentrations for Cl^- [113], and the unspecified ion X^- can be seen as essentially taking the role that Cl^- has in real systems.

Power spectrum analysis.

The power spectra (Fig 8) were computed with the fast Fourier-transform in MATLAB (http://se.mathworks.com/) and filtered to give one value per 0.1 log unit of the frequency.

Pyramidal cell model.

As neural model, we used the thick-tufted layer 5 pyramidal cell model by Hay et al. [62], which was implemented in the NEURON simulation environment [60]. The model was morphologically detailed (it had 196 sections, each of which we divided into 20 segments), and had a vertical extension of slightly less than 1300 μm from the tip of the basal dendrite to the tip of the apical dendrites. It contained ten active ion channels with different distributions over the somatodendritic membrane, including two Ca²⁺-channels (i^CaT, i^CaL), five K⁺-channels (i^KT, i^KP, i^SK, i^Kv3.1, i^M) and two Na²⁺-channels (i^NaT, i^NaS). In addition, it included a non-specific ion channel (i^h) and the non specific leakage current i^leak. The neuron had a membrane capacitance of 1 μF/cm² in the soma, and 2 μF/cm² in the dendrites, and leak conductances ranging between 0.0325 and 0.0589 mS/cm² over the somatodendritic membrane. We refer to the original publication for further model details [62].

Synapse model.

The neurons received Poissonian input trains through 10,000 synapses per neuron, a typical number for cortical neurons [114]. Each synapse had a mean input spike rate of 5 Hz. The synapses were uniformly distributed across the membrane so that the expected number of synapses in a segment was proportional to its membrane area. A population of ten neurons was simulated by running 10 independent simulations with the same neural model. The synapse distribution and spike trains were regenerated for each of the independent simulations (but with the same statistics in each case).

The synapses were modelled as α-shaped synaptic conductances: (24) where t₀ represents the time of onset. The time constant was set to τ = 2.0 ms. The maximal conductances of the synapses were set to g_max = 0.042 nS. For this value, the input evoked an average single-neuron action-potential (AP) firing rate of about 5 APs per second, which is a typical firing rate for cortical neurons [35, 63].

Population output to ECS.

The spatial extension of the cell morphology [62] was such that the maximal spatial distance between two segments (from tip of basal dendrite to tip of apical dendrite) was less than 1300 μm. We therefore considered a tissue depth of 1500 μm, and subdivided it into N = 15 ECS subvolumes (depth intervals) of length l_c = 100 μm, so that the neurons occupied the interior 13 subvolumes. Each neural segment was assigned as belonging to a particular subvolume n, determined by the spatial location of the segment midpoint. In the setup, the soma was placed in subvolume n = 3, the basal dendrites were in subvolumes n = 2, 3 and 4, and the apical dendrites were in subvolumes n = 3, …, 14. The multicompartmental model also included a short axon, which was, however, not based on the reconstruction and hence had no fixed coordinates. We assigned the axonal segments into the same subvolume as the soma, n = 3. The boundary subvolumes 1 and 15 contained no neural segments (see Fig 1A).

The transmembrane current density () of ion species k is available in the NEURON simulation environment. It was multiplied by the surface area of the segment (A_seg) to get the net current, and divided by Faraday’s constant (F) to get a net ion flux with units mol/s. During the neural simulation, we grouped all currents that were carried by a specific ion species into the net transmembrane influx/efflux of this ion species. We assumed that all non-specific currents, including the synaptic currents (i^leak, i^h, i^syn) were carried by a non-specific anion that we denoted X^-. In this way we could compute the net efflux of each ion species into a subvolume n: (25) Here, the sum was taken over all neural segments (seg) of all 10 neurons contained in subvolume n. The factor 2 in the denominator in the expression for was due to Ca²⁺ having valence 2, and the negative sign in the expression for was due to X⁻ having valence -1. We also kept track of the (non-ionic) capacitive currents, as required by the electrodiffusive formalism (Eq 14): (26) The intracellular dynamics was directly adopted from the original model [62]. Transmembrane currents there had no effect on intracellular ion concentrations, except for Ca²⁺ concentration, which was modelled to account for Ca²⁺ dependent K⁺ channels.

For technical reasons (concerning memory usage), 84 s of output from a single neuron was generated in the following way: First, we ran ten simulations, each producing 10 s of activity. Next, we removed the initial 1.6 s of all the ten simulations, to remove the transient neuronal activity observed initially in the simulations, leaving us with ten 8.4 s time series. Finally, these were used as successive output periods from a single neuron, and combined into a total 84 s time series.