To calculate the extracellular electric potential (LFP) we will make use of a recently developed kernel method based on unitary-LFP's (uLFP) (Telenczuk et al., 2020). This method was designed for the estimation of LFP's from spiking network simulations and consists in convolving the spikes of the network with two kernels, according to the formula:

where Ke(x→,t-τ) and Ki(x→,t-τ) are the kernels associated with spikes from excitatory and inhibitory neurons, respectively, while {t_e,j} and {t_i,j} are the spiking times of excitatory and inhibitory neurons. In mean-field models, we typically deal with the mean firing activity of two populations, excitatory and inhibitory neurons, which we will note by ν_e and ν_i. To adapt this method to mean-field models, the precise timing spikes in the previous expression can be replaced by the rates of spiking activity, ν_e and ν_i, which yields:

where we have assumed that the spiking rates are high enough such that the typical inter-spike interval is small compared to the characteristic time of the kernel and thus the discrete spiking times can be replaced by the average firing rate. We notice that a typical size of a neuronal population in a mean-field is of thousands of cells with individual firing rates up to tens of Hz. Thus, for an electrode registering the activity of the population, this gives an inter-spike interval in the order of 1x10⁻² ms, while the kernels (as shown below) have characteristic times of a few ms, which validates the previous assumption.

The kernels functions K_e and K_i were estimated from experimental data, by fitting a Gaussian template to the relation between single spikes and LFP's obtained in micro-electrode recordings (i.e., unitary LFP's) (Teleńczuk et al., 2017; Telenczuk et al., 2020). The kernels at position x→ and time t are given by:

where A is the amplitude (which can be negative), σ is the standard deviation of the kernel in time, and t_p is the peak time of the kernel. The latter is given by

where t₀ is the time of the spike of the cell, |x→-x→0| is the distance between position of the electrode x→ and the position x→0 of the current source, d is a constant delay, and v_a is the axonal speed. We use the value of v_a = 200 mm/s, estimated from human LFP recordings (Teleńczuk et al., 2017).

To model the observed near-exponential amplitude decay with distance, the following expression can be used for A(x→), in cylindrical coordinates (Teleńczuk et al., 2017):

where A₀(z) is the maximal amplitude, which depends on coordinate z (cortical depth), |x − x₀| is the radial distance between position of the electrode x→ and the source at x→0 (see coordinate system in Figure 1 for a cylindrical region) and λ is the space constant of the decay. From human microelectrode recordings, λ was consistently found around 200–250 μm (Teleńczuk et al., 2017). A diagram of the kernel and the corresponding values for A₀(z) are shown in Figure 1.

To perform the calculation of the LFP from mean-field models we will replace the single neuron magnitudes in the previous expressions by their mean values over the neuronal population. For the spatial dependence we take the average value of the function A(x→) over a circular region around the point x→ where the field is measured. Taking into account the exponential decay, we will assume that the radius of the region of interest around x→ (r_region) is defined by the characteristic distance λ, and we will take r_region = 2λ. For the temporal part of the kernel we will neglect the term corresponding to the axonal propagation, which holds valid as far as the propagation time is small compared with the typical temporal scale of the kernel, for example in the case where the region comprised by the neuronal population is small enough, as we will show in the example provided in the Results section. Thus, we will take: A(x→)≈A0(z)<exp[-|x-x0|/λ]>|r=2λ=A0(z)12(1-3e2)≈0.3A0(z) (independent of λ) and t_p = d, where ‘<>’ indicates the mean over the area defined by the radius r_region = 2λ.

To estimate the magnetic field generated by a neuronal population we will assume that the variations in charge and current densities are slow enough (frequencies below the thousands of hertz) such that time derivatives in Maxwell equations can be neglected. Within this quasi-static approximation, the magnetic field generated by a current density J(r, t) is given by (Geselowitz, 1970; Hagen et al., 2018; Garcia-Rodriguez and Destexhe, 2021):

where μ₀ is the magnetic susceptibility, r′ is the location of the source and r the point where the field is measured. For the case where |r-ro′|>>|r′-ro′|, where ro′ represents the location of the center of the neuronal population we can take the far-field approximation and we have:

where Q = ∫J(r′, t)d³r′ is the total dipole moment. This is in particular valid for the case of MEG, where sensors are located outside the scalp, far from the neuronal population generating the field.

In the brain, neurons are surrounded by the extracellular medium which is an electrical conductor. The current J is then composed of three main sources: synaptic currents (J^s), action potentials (J^a) and conduction currents (J^c), the last ones being related to the induced currents in the conducting surrounding medium. Thus, we can write Q(r, t) = Q^s + Q^a + Q^c. Currents related with action potentials tend to cancel each other, and their contribution to the magnetic field can be neglected (Lindén et al., 2014). In the following we show how to estimate the contributions from the conduction and synaptic currents. The conduction currents, given a medium conductivity σ(r, t), can be written as:

where E(r, t) is the electric field at point r and Φ is the local electric field potential (LFP). The conductivity σ(r, t) will in general depend on the position and may vary with time. We will compartmentalize the brain or region under study into small volume domains of constant conductivity, such that we can write for Q^c (Geselowitz, 1970; De Munck et al., 2012):

where the summation goes over all surfaces of discontinuity, σi(t)′ and σi(t)″ are the conductivities at the two regions separated by the surface i, ds′ is the differential area element for surface i and where the identity ∫Vi∇′Φ(r′,t)d3r′=∫S(i)Φ(r′,t)dS′i was used. We notice that the electric field potential Φ can be estimated from the mean-field models using the kernel method described in the previous section, from where the contribution of Q^c to the magnetic field is fully determined. We notice in addition that the conductivity σ may also depend on the frequency and direction. For simplicity in this paper we will assume an isotropic medium and independent of the frequency, but in principle these dependencies can be incorporated within our method.

To estimate the contributions of intrinsic neuronal currents to the MEG we will take into account that axial currents are believed to be the main source for the magnetic-field (in opposition to transmembrane currents) (Hämäläinen et al., 1993). We notice however that, in general, the mean-field and AdEx models are based on point-neurons (i.e., a single-compartment neuronal model). In order to incorporate axial currents to our model we will extend our description and consider the mean-field model in combination with a two compartment model (sometimes referred to as an “hybrid modeling”; Hagen et al., 2018), which is the minimum configuration that incorporates intra-neuron current flow. We show the diagram of the two-compartment model in Figure 1A. One of the compartments corresponds to the soma-perisomatic region and the other to the apical dendrites. We will assume that the membrane voltage for each of these compartment is given by:

where g_L and E_L are the conductance and reversal potential of the leakage channel, Isynj is the synaptic input to compartment j, W is a variable describing neuronal adaptation (which we only consider affecting the first compartment for simplicity) and I_A is the axial current between the two compartments, i.e., I_A = (V₁ − V₂)/R_A, being R_A the axial resistance.

Following Kuhn et al. (2004), we will assume that the mean membrane potential of each compartment can be obtained by taking the stationary solution of (Equation 10) to static synaptic currents given by the synaptic bombardment with firing rates ν_e and ν_i (Di Volo et al., 2019). Thus, subtracting both expressions in Equation 10, and taking the stationary solution we get (for g_L1 = g_L2 = g_L):

where μ_Vj indicates the mean voltage of compartment j and from where the mean axial current (over the neuronal population) can be estimated. We notice that in general the synaptic currents will be a function of the voltage of each compartment and the firing rates of the population, in which case both stationary solutions of Equation 10 have to be considered. In a very general way we can write Isynj=Isynexj+Isyninj, with Isynexj,Isyninj the excitatotory and inhibitory synaptic currents, respectively, given by:

where μ_Ge(i) is the mean excitatory (inhibitory) synaptic conductance and Ke(i)j is the number of excitatory (inhibitory) synaptic connections arriving to compartment j (Di Volo et al., 2019). The estimation of the mean synaptic conductances is to be defined from the mean-field model of use. In the examples presented in this paper we will consider a conductance based network, where the conductance at each synaptic terminal is increased by a quantal amount q_e(i) for each pre-synaptic spike (see Section 3 for details). For these networks and assuming that the spiking activity follows a poissonian distribution, then the mean conductances can be written as (Di Volo et al., 2019):

where τ_e(i) is the characteristic decay time of the synaptic conductance. With the definition of the synaptic currents and mean conductances, Equations (10) and (11) are fully determined and the axial currents can be calculated (see Supplementary material for further details). By further assuming that excitatory neurons are the main source of the magnetic fields then we have:

which gives the mean dipole moment generated by intrinsic neuronal currents, where L is a characteristic length of the dipole and n_e is the number of excitatory neurons. Thus, together with Equation (9) (Q^c) and (7) the magnetic field is completely determined.

In the Results section we will present an example of the use of the method for a recently developed mean-field model of a network of Adaptive-exponential-Integrate-and-Fire (AdEx) neurons.

The Virtual Brain simulator platform (TVB) provides an environment for large-scale brain simulations, where the brain is compartmentalized in macro-domains represented as nodes in a network. Each node in the network is described via a mean-field model. The connectivity among brain regions is provided via a brain connectome. The size of the compartmentalization can be defined by the user, depending on the characteristic of the study and computational resources. For our simulations we will consider a coarse-grain simulation where each node represents a different brain region. For more details on our simulations see the Section 3. For details on the TVB platform see Sanz Leon et al. (2013) and Melozzi et al. (2017).

To provide an example and validate our method we will consider a recently developed mean-field model of a network of inhibitory and excitatory Adaptive-Exponential-Integrate-and-Fire (AdEx) neurons (Brette and Gerstner, 2005). The spiking AdEx neuronal model is defined by the system of equations:

where V is the membrane potential, w is an adaptation current, C = 200 pF is the membrane capacity, g_L = 10 nS is the leakage conductance, E_L = −63 mV is the leakage reversal potential, V_T = −50 mV, Δ = 2 mV (0.5 mV) for excitatory (inhibitory) neurons, I_syn is the synaptic current, a is the sub-threshold adaptation constant and b is the spiking adaptation constant. When V > V_T at time t = t_sp a spike is generated, the membrane potential is reset to V_res = −65 mV and remains at this value for a refractory time t_ref = 5 ms, and the adaptation variable is increased by an amount b. In our simulations we will consider that no adaptation occurs in inhibitory neurons (i.e., a = b = 0), while for excitatory neurons we will consider b = 60 pA and a = 0. The synaptic current I_syn received by a neuron is given by the spiking activity of all presynaptic neurons connected to it. The total synaptic current can be written as a sum of the excitatory and inhibitory synaptic activity I_syn = G_e(E_e − V) + G_i(E_i − V), where E_e = 0 (E_i = −80 mV) is the excitatory (inhibitory) reversal potential and G_e, G_i are the synaptic conductances. We model the synaptic conductances as decaying exponential functions that experience a quantal increase q_e and q_i at each pre-synaptic spike: Ge(i)=qe(i)∑Θ(t-tsp)et-tspτe(i), where q_e = 1.5 nS, q_i = 5 nS, τ_e = τ_i = 5 ms.

We will consider a network of 10,000 neurons, with 80% of excitatory and 20% of inhibitory neurons. Neurons in the network are randomly connected with probability p = 5%.

The mean-field equations for the AdEx network are given to a first-order by (Di Volo et al., 2019):

where ν_e,i is the mean neuronal firing rate of the excitatory and inhibitory population, respectively, W is the mean value of the adaptation variable, F is the neuron transfer function (i.e., output firing rate of a neuron when receiving excitatory and inhibitory inputs with mean rates ν_e and ν_i and with a level of adaptation W), a and b are the sub-threshold and spiking adaptation constants, t_w is the characteristic time of the adaptation variable, T is a characteristic time for neuronal response, μ_V is the average membrane voltage and E_L is the leakage reversal potential (see Di Volo et al., 2019 for details).

To validate the results of our method for the calculation of magnetic fields, we will compare the estimation from our method with the results obtained from a two-compartment spiking AdEx network. In this network inhibitory neurons will be consider as point neurons as before, while excitatory neurons will be described by a two-compartment model in analogy to the model proposed in Equation (10). The equations corresponding to the peri-somatic area (Equations 18 and 19) will take into account the spiking mechanism and the adaptation. We will consider that the soma receives inputs from synaptic currents and from the axial current coming from the dendritic region:

Where C_S = 200 pF, g_LS = 10 nS, E_L = −63 mV, V_T = −50 mV, Δ_T = 2 mV, τ_w = 500 ms, b = 60 pA a = 0.

The axial current is defined as before I_A = g_A(V_D−V_S) with g_A = 400nS.

The dendritic membrane potential is described by a leak current and receives input from synaptic currents. The membrane voltage of the dendrite is described by:

With C_D = 10 pF, and g_L−D = 2 nS.

For both populations when an action potential is emitted (voltage is greater than V_T), the system is reset to its resting value and remains at that value for a refractory period of t_ref = 5 ms. We will take V_T = −40 mV, V_R = −55 mV for excitatory cells and V_T = −47.5 mV, V_R = −65 mV for inhibitory cells.

To estimate the magnetic field from the spiking network we will take the mean axial current generated in the excitatory neurons of the network and proceed to calculate the mean dipole moment as in Equation (14).

We provide next two examples of the method to calculate the LFP. First we show the results for a single mean-field model of an AdEx network, and the second example consists in a large-scale simulation of the human brain performed with The Virtual Brain simulator platform.

In Figure 2 we show the results of the mean-field model and the LFP calculations for the first example. To validate our results we show the comparison with the calculation of the LFP for the corresponding spiking AdEx network. For the latter case we show the raster plot and the corresponding mean firing rate of the excitatory and inhibitory populations. We can see that the results from the mean-field model can correctly capture the amplitude and duration of the up-and-down states obtained in the LFP from spiking-networks. The average amplitude of the high-state in the LFP is of (206 ± 26) μV for the spiking case and (209 ± 2) μV for the mean-field case. The average duration of the high-state in the LFP is (0.59 ± 0.14) s for the method applied to the spiking network and (0.57 ± 0.09) s for the mean-field case. The duration is measured as the time between the beginning and the end of the high state, where the criteria for the starting and ending points is defined as the crossing point of 3 standard deviations from the base state.

The second example we present consists in a large-scale simulation of the human brain, performed with the TVB simulator platform. In this simulation each region of the brain is represented by one neuronal ensemble of AdEx neurons (as described in the first example). The connection between regions was defined by human tractography data (https://zenodo.org/record/4263723,Berlinsubjects/QL_20120814) from the Berlin empirical data processing pipeline (Schirner et al., 2015; Goldman et al., 2021). A parcellation of 68 regions is used, with long-range excitatory connections and delays defined by 120 tract length and weight estimates in human diffusion magnetic resonance imaging (dMRI) data (Sanz-Leon et al., 2015; Goldman et al., 2021). The LFP is calculated using the kernel-method, for which it has been incorporated in the TVB platform. For details on the connectome and the simulations with the TVB see Sanz Leon et al. (2013), Schirner et al. (2015), Sanz-Leon et al. (2015), and Goldman et al. (2021). The results are shown in Figure 3, where we show only the LFP for the Surface depth. All the simulations presented in the paper were performed using the python language and the TVB simulation platform for whole brain simulations.

In the following we proceed with the calculation of the magnetic field for the example of the AdEx mean-field model, as explained in the Section 2. As for the LFP, we consider two examples, one consisting in a single mean-field and another comprising a large-scale simulation of the human brain.

To perform the calculations two points are to be specified: i) the synaptic input and ii) the specific geometry and neuronal distribution of the area under study. The first one is defined by the mean-field model adopted for the simulations while the second depends on the specific features of the study. In our case the synaptic currents are given by Equation (12) and (13). For our simulations we adopt τ_e(i) = 5ms, q_e(i) = 1.5nS (5 nS) and E_e(i) = 0 (−80 mv). For Ke(i)j we take Ke(i)j=ne(i)pe(i)j with n_e(i) the number of excitatory (inhibitory) neurons in the population and pe(i)j the connectivity probability of an excitatory (inhibitory) synapse at compartment j. To define pe(i)j we assume as before that every two neurons in the network are connected with a probability p_o = 0.05. In addition, following experimental data, we will assume that inhibitory synapses have a higher concentration around the soma while excitatory synapses dominate in the apical dendrites (Megıas et al., 2001). We will consider a 60–40% distribution for inhibitory synapses (60% in the peri-somatic region, 40% at the apical dendrites) and a 30–70% distribution for excitatory synapses (see Figure 4A for a diagram of the two-compartment model). Thus, pe(i)1=0.3po (0.6p_o) and pe(i)2=0.7po (0.4p_o). Given the synaptic currents for each compartment, the axial current and its associated magnetic field is determined from Equation (7), (10), and (14).

To perform the simulations a concrete geometry and neuronal distribution for the region under study must be defined. For this, we will assume that neurons are aligned as shown in Figure 4B and uniformly distributed. It is easy to see that current dipoles oriented perpendicular to the scalp have a smaller or null contribution to the magnetic field and that the MEG is rather driven by dipoles oriented tangentially, which is a commonly accepted feature (Hämäläinen et al., 1993). For our simulations we will consider cubic regions as indicated in Figure 4C and for simplicity we will assume that only neurons and dipole moments with purely tangential currents (parallel to the scalp) contribute to the magnetic field. Furthermore, we will assume that the vector r-ro′ between the dipole and the sensor is aligned in the radial direction (i.e., perpendicular to the scalp), such that the calculated field corresponds to a purely tangential component (a generalization for an arbitrary angle is straightforward). We assume that the sensor is located at a distance of 3 cm from the neuronal population. For the contribution of the conduction currents (Equations 8 and 9) we will consider a single interface between two regions of different conductivity. In particular we will assume that the two conductivities correspond to the one of the gray matter and the cerebrospinal fluid, respectively, which exhibit a relatively high difference in conductivity (we use σ = 0.3 and 2.1 S/m, respectively, McCann et al., 2019; Mandija et al., 2021). The calculation of the LFP for Equation (9) and the related parameters follow the ones indicated for Figure 2.

In Figure 5B we show the results of our simulations. We show here the total magnetic field and the field generated by each different source (synaptic and conductance currents). For the conductance currents we show the contribution of each depth to the field. To validate our results we show the results obtained from the spiking AdEx two-compartment model in Figure 5A. For the spiking networks we show the raster plot, the mean firing rates and the resulting magnetic field. For this last one we only consider the contribution of synaptic (i.e., axial) currents to the field. The results from the mean-field method can reasonably capture the amplitude and duration of the up-and-down states obtained from the spiking network. The average amplitude of the up-state in the MEG is of (276 +/− 40) fT for the spiking case and (305 +/− 3) fT for the mean-field case. The average duration is of (0.87 +/− 0.43) s for the spiking case and (0.53 +/− 0.01) s for the mean-field case.

From these simulations we can see that the axial neuronal currents constitute the main source for the MEG, which is in agreement with typical experiments. In our simulations we took only a single interface between regions of different conductivity. The contribution of the different sources to the MEG will depends of course on the particular geometry, conductivities and neuronal distributions. Nevertheless, it is usually accepted that the magnetic field is mainly driven by axial neuronal currents (Hämäläinen et al., 1993). Our analysis would agree with this, but it also provides means to further analyze the different contributions and in particular to explore the possible contributions of different layers. For a comparison with experimentally measured slow-waves with MEG see for example Simon et al. (2000).

To end this section we present the results of the MEG calculations for a large-scale simulation of the human brain with The Virtual Brain platform, corresponding to the same set of simulations as in Figure 3. For the MEG we aim to describe the fields measured by N sensors and generated by the M sources (brain regions in our case). The magnetic field measured by a sensor at position r outside the conductor can be written in general as (Mosher et al., 1999 and Sanz-Leon et al., 2015):

where Q indicates the corresponding dipole moments, r_Q is the location of the sources (dipoles) and G(r, r_Q) is called the gain or transfer matrix and contains all the information regarding the geometry of the specific head-model under consideration. This is a generalization of Equation (7). In the TVB platform the gain matrix can be computed and registered for different head-models given the locations of sources and sensors (Sanz-Leon et al., 2015). In our simulations we will consider that each dipole moment is located at the center of the corresponding brain region and that there is one sensor per each brain region. For simplicity we will consider G=μ04π|r−r′|2I, where, as before, we assume that the vector r − r′ is aligned in the radial direction, |r − r′| is a typical distance between the sensor and the source, and I is the identity matrix, such that each sensor measures only the field generated by each brain region. For a more detailed calculation of MEG, a more detailed head-model and orientation of dipoles should be taken into consideration.

We show in Figure 6 the results of our whole brain simulations. In Figure 6A of the figure we show the time series of the mean firing rates of each region and the corresponding magnetic field for each sensor registered in the simulation. For these simulations we are assuming an homogeneous medium (i.e., uniform conductivity), from where the only source of the magnetic field are the neuronal axial currents. In Figure 6B we show a colormap of the magnetic field for region corresponding to a single slow-wave in our simulations.