Skip to main content

METHODS article

Front. Comput. Neurosci., 13 October 2022

Mean-field based framework for forward modeling of LFP and MEG signals

  • CNRS, Paris-Saclay Institute of Neuroscience (NeuroPSI), Paris-Saclay University, Saclay, France

The use of mean-field models to describe the activity of large neuronal populations has become a very powerful tool for large-scale or whole brain simulations. However, the calculation of brain signals from mean-field models, such as the electric and magnetic fields, is still under development. Thus, the emergence of new methods for an accurate and efficient calculation of such brain signals is currently of great relevance. In this paper we propose a novel method to calculate the local field potentials (LFP) and magnetic fields from mean-field models. The calculation of LFP is done via a kernel method based on unitary LFP's (the LFP generated by a single axon) that was recently introduced for spiking-networks simulations and that we adapt here for mean-field models. The calculation of the magnetic field is based on current-dipole and volume-conductor models, where the secondary currents (due to the conducting extracellular medium) are estimated using the LFP calculated via the kernel method and the effects of medium-inhomogeneities are incorporated. We provide an example of the application of our method for the calculation of LFP and MEG under slow-waves of neuronal activity generated by a mean-field model of a network of Adaptive-Exponential Integrate-and-Fire (AdEx) neurons. We validate our method via comparison with results obtained from the corresponding spiking neuronal networks. Finally we provide an example of our method for whole brain simulations performed with The Virtual Brain (TVB), a recently developed tool for large scale simulations of the brain. Our method provides an efficient way of calculating electric and magnetic fields from mean-field models. This method exhibits a great potential for its application in large-scale or whole-brain simulations, where calculations via detailed biological models are not feasible.

1. Introduction

The electric and magnetic fields generated by neuronal currents are two commonly used brain signals to study neuronal activity. The measurement of these two signals involve different experimental techniques that can capture neuronal activity at different scales. The registration of electric local field potentials (LFP) involves invasive experimental techniques to measure the electric potential in the cerebral extracellular medium and it captures the activity of thousands of nearby neurons, in particular their synaptic activity (Niedermeyer and Lopes Da Silva, 2005; Buzsáki et al., 2012). On the other hand, the measurement of the magnetic field in MEG (magnetoencephalogram) is a non-invasive technique that captures the field generated by thousands to millions of neurons. Among non-invasive techniques to measure neuronal activity, MEG provides a relatively high spatial and temporal resolution, which are currently in the order of millimeters and milliseconds, respectively (De Pasquale et al., 2010; Hansen et al., 2010; Stokes et al., 2015).

The calculation of these signals from the underlying neuronal sources (usually referred to as “forward-modeling”) can be performed via detailed biophysical models (Hämäläinen et al., 1993; Bédard et al., 2004; Nunez et al., 2006; Bédard and Destexhe, 2012; Lindén et al., 2014; Hagen et al., 2018; Ilmoniemi and Sarvas, 2019). Although successful, this turns out to be very computationally demanding, specially when large populations of neurons are into consideration. Nowadays the utilization of mean-field models describing the activity of neuronal populations has become a powerful tool for large-scale simulations (Sanz Leon et al., 2013; Depannemaecker et al., 2021), however the calculation of brain signals from the mean-field models is still under development. Thus, the emergence of new modeling-frameworks that allow an accurate and efficient calculation of multiple brain-signals from mean-field models is currently of great interest.

In this paper we provide a method to calculate the LFP and MEG signals from mean-field models. To calculate the LFP we base our method on the unitary LFP (uLFP), which is the LFP generated by the post-synaptic currents from a single axon. We adopt a recently developed kernel-method based on uLFP's for spiking neural networks (Telenczuk et al., 2020) and we adapt this method for its application to mean-field models. Kernel methods based on neuronal models (Hagen et al., 2016) and experimental data (Telenczuk et al., 2020) have been proposed for the calculation of LFP's. In our approach we use a kernel method based on experimental data of uLFP's (Telenczuk et al., 2020), which is not restricted by specific model constraints.

To calculate the magnetic-field we make use of current-dipole and volume-conductor models for an inhomogeneous medium. We estimate the contribution of neuronal currents to the magnetic field from a mean-field model and we incorporate the previous estimation of the LFP to calculate the contribution of secondary currents due to the medium conduction. We provide an example of our method for the calculation of the LFP and MEG under slow-wave activity generated by a mean-field model of Adaptive-Exponential Integrate-and-Fire (AdEx) neurons. To validate our method we compare the results with the ones obtained from simulations of the corresponding spiking neural networks. Finally we provide an example of our method for whole brain simulations performed with The Virtual Brain (TVB), a recently developed tool for large scale simulations of the brain (Sanz Leon et al., 2013; Goldman et al., 2020).

2. Method

2.1. Calculation of the electric potential from mean-field models

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:

Ve(x,t)=  Ke(x,tτ)(jδ(τte,j))dτ                    +  Ki(x,tτ)(jδ(τti,j))dτ ,    (1)

where Ke(x,t-τ) and Ki(x,t-τ) are the kernels associated with spikes from excitatory and inhibitory neurons, respectively, while {te,j} and {ti,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:

Ve(x,t)= Ke(x,tτ) νe(τ) dτ                    + Ki(x,tτ) νi(τ) dτ.    (2)

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−2 ms, while the kernels (as shown below) have characteristic times of a few ms, which validates the previous assumption.

The kernels functions Ke and Ki 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:

K(x,t)=A(x) exp[-(t-tp)2/(2σ2)],    (3)

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

tp=t0+d+|x-x0|/va,    (4)

where t0 is the time of the spike of the cell, |x-x0| is the distance between position of the electrode x and the position x0 of the current source, d is a constant delay, and va is the axonal speed. We use the value of va = 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):

A(x)=A0(z) exp[-|x-x0|/λ],    (5)

where A0(z) is the maximal amplitude, which depends on coordinate z (cortical depth), |xx0| is the radial distance between position of the electrode x and the source at x0 (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 A0(z) are shown in Figure 1.

FIGURE 1
www.frontiersin.org

Figure 1. (A) Diagram of a pyramidal neuron and the coordinate system. The dashed lines indicate the two compartments considered for the calculation of the magnetic fields. (B) Experimental measurements of the uLFP (circles) at different distances from the soma (x-coordinate) together with the uLFP kernel of Equation (3) (solid line). (C) Numerically estimated uLFP amplitude as a function of depth (z-coordinate) for inhibitory neurons (white and black circles correspond to different simulations of the same system, see Telenczuk et al., 2020). (D) Amplitude of the uLFP for inhibitory and excitatory neurons at selected depths: Deep layers (z = –400 μm), Soma (z = 0),Superficial layers (z = 400 μm) and Surface (z = 800 μm). Adapted from Telenczuk et al. (2020).

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 (rregion) is defined by the characteristic distance λ, and we will take rregion = 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 tp = d, where ‘<>’ indicates the mean over the area defined by the radius rregion = 2λ.

2.2. Calculation of the magnetic field from mean-field models

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):

B(r,t)μ04π J(r,t)×rr|rr|3d3r    (6)

where μ0 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:

B(r,t)μ04πQ×rro|rro|3.    (7)

where Q = ∫J(r′, t)d3r′ 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 (Js), action potentials (Ja) and conduction currents (Jc), the last ones being related to the induced currents in the conducting surrounding medium. Thus, we can write Q(r, t) = Qs + Qa + Qc. 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:

Jc=σ(r,t)E(r,t)=-σ(r,t)Φ    (8)

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 Qc (Geselowitz, 1970; De Munck et al., 2012):

Qc=σ(r,t)Φ(r,t)d3r       =i(σi(t)σi(t))Φ(r,t)dSi    (9)

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)dSi 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 Qc 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:

C1dV1dt=gL1(ELV1)+Isyn1W+IAC2dV2dt=gL2(ELV2)+Isyn2IA    (10)

where gL and EL 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 IA is the axial current between the two compartments, i.e., IA = (V1V2)/RA, being RA 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 gL1 = gL2 = gL):

μV1-μV2=Isyn1-Isyn2-W2/RA+gL    (11)

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:

Isynexj=KejμGe(Ee-μVj)Isyninj=KijμGi(Ei-μVj)    (12)

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 qe(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):

μGe(i)=νe(i)τe(i)qe(i)    (13)

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:

Qs=neLIAe    (14)

which gives the mean dipole moment generated by intrinsic neuronal currents, where L is a characteristic length of the dipole and ne is the number of excitatory neurons. Thus, together with Equation (9) (Qc) 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.

2.3. The TVB simulation platform

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).

3. Results

We present in the following an example of our method for the calculation of the Local Field Potential (LFP) and magnetic fields from a mean-field model. We will validate our method via comparison with results obtained from simulations of the corresponding spiking neural networks. The example we present is based on simulations of slow-waves of neuronal activity obtained with a mean-field model of Adaptive-Exponential Integrate-and-Fire neurons. Slow-waves (0.5-3.5 Hz, δ-band) are a key feature of neuronal activity during deep-sleep and can be experimental measured via LFP and MEG (Destexhe et al., 1999; Simon et al., 2000; Lopes da Silva, 2013).

3.1. AdEx mean-field model

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:

CdVdt=gL(EL-V)+gLΔ exp(V-VTΔT)-w+Isyn    (15)
dwdt=bδ(t-tsp)+aτw(V-EL)-wτw    (16)

where V is the membrane potential, w is an adaptation current, C = 200 pF is the membrane capacity, gL = 10 nS is the leakage conductance, EL = −63 mV is the leakage reversal potential, VT = −50 mV, Δ = 2 mV (0.5 mV) for excitatory (inhibitory) neurons, Isyn is the synaptic current, a is the sub-threshold adaptation constant and b is the spiking adaptation constant. When V > VT at time t = tsp a spike is generated, the membrane potential is reset to Vres = −65 mV and remains at this value for a refractory time tref = 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 Isyn 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 Isyn = Ge(EeV) + Gi(EiV), where Ee = 0 (Ei = −80 mV) is the excitatory (inhibitory) reversal potential and Ge, Gi are the synaptic conductances. We model the synaptic conductances as decaying exponential functions that experience a quantal increase qe and qi at each pre-synaptic spike: Ge(i)=qe(i)Θ(t-tsp)et-tspτe(i), where qe = 1.5 nS, qi = 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):

Tdνe,idt=Fe,i(W,ν¯e,νi)νe,idWdt=Wτw+bνe+a(μV(ν¯e,νi,W)EL)    (17)

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, tw is the characteristic time of the adaptation variable, T is a characteristic time for neuronal response, μV is the average membrane voltage and EL is the leakage reversal potential (see Di Volo et al., 2019 for details).

3.2. Two compartment AdEx model

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:

CSdVSdt=gLS(EL-V)+gLSΔTexp(V-VTΔT)-w+IA+Isyn-S    (18)
dwdt=bδ(t-tsp)+aτw(VS-EL)-wtw    (19)

Where CS = 200 pF, gLS = 10 nS, EL = −63 mV, VT = −50 mV, ΔT = 2 mV, τw = 500 ms, b = 60 pA a = 0.

The axial current is defined as before IA = gA(VDVS) with gA = 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:

CDdVDdt=gLD(EL-VD)-IA+Isyn-D    (20)

With CD = 10 pF, and gLD = 2 nS.

For both populations when an action potential is emitted (voltage is greater than VT), the system is reset to its resting value and remains at that value for a refractory period of tref = 5 ms. We will take VT = −40 mV, VR = −55 mV for excitatory cells and VT = −47.5 mV, VR = −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).

3.3. LFP from the AdEX mean-field model

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.

FIGURE 2
www.frontiersin.org

Figure 2. (A) LFP calculations obtained with the kernel method on spiking-network simulations (in blue) at different depths (Surface, Superficial, Soma and Deep, see table in Figure 1). Simulations were performed on a network of 10,000 AdEx neurons and LFP was calculated according to Equation (1). Neurons where uniformly distributed over a square surface of 1 x 1 mm2 (see diagram in Figure 4C for reference). We show on the top the raster plot of the neuronal activity during the simulation, for excitatory (green) and inhibitory (red) neurons together with the corresponding firing rate for each population. (B) LFP calculations obtained with the kernel method from mean-field simulations. We show the average firing rate for the excitatory and inhibitory neurons (top) in a network of 10,000 AdEx neurons, computed from Equation (17) and the corresponding LFP obtained with Equation (2) (blue).

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.

FIGURE 3
www.frontiersin.org

Figure 3. Whole brain simulations of the human brain performed with The Virtual Brain (TVB) platform. (A) Top: time course of the average firing rates for each brain region in the simulations. Bottom: Local field potential calculated with the kernel method for each region. (B) Colormap of the LFP amplitude at each brain region for a down (t = 0) and up (t = 0.26 s) states during one slow-wave oscillation. For details on the TVB simulator of the human brain see Sanz Leon et al. (2013), Schirner et al. (2015), Sanz-Leon et al. (2015), and Goldman et al. (2021).

3.4. Magnetic field from the AdEx mean-field model

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, qe(i) = 1.5nS (5 nS) and Ee(i) = 0 (−80 mv). For Ke(i)j we take Ke(i)j=ne(i)pe(i)j with ne(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 po = 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.6po) and pe(i)2=0.7po (0.4po). Given the synaptic currents for each compartment, the axial current and its associated magnetic field is determined from Equation (7), (10), and (14).

FIGURE 4
www.frontiersin.org

Figure 4. (A) Diagram of the two compartment neuronal model used for MEG calculations. One compartment covers the soma and perisomatic region with high concentration of inhibitory synapses (red) and the second compartment covers the apical dendrites with dominant concentration of excitatory synapses (green). (B) Diagram of the cerebral cortex with the arrangement and orientation of neurons in the gyrus and sulcus. The different depth layers used for LFP calculations are indicated by dashed lines. (C) Diagram of the domains, layers and surface consider for the MEG simulations. Xm indicates the position of the sensor.

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.

FIGURE 5
www.frontiersin.org

Figure 5. Calculation of the MEG signal. (A) Simulations from a spiking two-compartment AdEx network. We show the raster plot (top), the mean firing rate (middle) and the magnetic field generated by axional neuronal currents (bottom). (B) Results obtained with our method from the AdEx mean-field model. In the first plot (top) we show the total magnetic field generated by the system (axial neuronal currents and conduction currents, dashed black line) together with the contribution of the axial currents (solid red line). Bottom: Contribution to the magnetic field by of the conduction currents. The contribution of the different depths is shown.

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):

B=G(r,rQ)Q    (21)

where Q indicates the corresponding dipole moments, rQ is the location of the sources (dipoles) and G(r, rQ) 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π|rr|2I, where, as before, we assume that the vector rr′ is aligned in the radial direction, |rr′| 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.

FIGURE 6
www.frontiersin.org

Figure 6. Whole brain simulations of the MEG signal in the human brain performed with The Virtual Brain (TVB) platform. (A) Top: time course of the mean firing rates for each brain region in the simulations. Bottom: magnetic field for each region calculated with the method introduced in this paper. (B) Colormap of the magnetic field amplitude at each brain region for a down (t = 0) and up (t = 0.26 s) states during one slow-wave oscillation. For details on the TVB simulations of the human brain see Sanz Leon et al. (2013), Schirner et al. (2015), Sanz-Leon et al. (2015), and Goldman et al. (2021).

4. Discussion

In this paper we have introduced a method for forward modeling of local field potentials (LFP) and magnetic fields generated by neuronal activity. Our method allows to calculate these signals from mean-field models of neuronal populations. We have proposed a method based on a phenomenologically defined kernel to calculate the LFP, which provides an estimation of the field as a function of depth and distances from the cell body. To calculate the magnetic-field we utilized a current-dipole and volume conductor model which we combined with the kernel method of LFP to estimate the secondary currents induced in the conducting medium by the neuronal electric fields. Our method considers a very general scenario with a non-homogeneous conducting medium. Other efforts to calculate LFP signals from kernel methods have been previously proposed (Hagen et al., 2016; Skaar et al., 2020), where kernels are estimated from simulations of detailed biophysical neuronal models. Our method is based on experimentally measured uLFP's, which is not limited by specific model constraints.

We have presented an example of the application of our method for a mean-field of a population of Adaptive Exponential Integrate and Fire neurons. For the example shown in the paper the magnetic-field is mainly generated by neuronal currents (i.e., primary currents) with a nearly negligible contribution of the conduction currents (secondary currents).

One possible source of limitations of our method is associated with the characteristic of the mean-field models. In particular, the mean-field models presented here correspond to randomly connected networks, which prevents the utilization of specific network structures.

Finally, our method exhibits a great potential for application in large-scale or whole-brain simulations, where calculations via detailed biological models are not feasible. Further improvement of our method would include the refinement of the kernel used for the LFP calculations, which is envisioned for the near future. The specification of the kernel for different brain-regions and a more detailed analysis of its depth-dependence will lead to more accurate calculations of the local field potential and the MEG signal.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. The code for the implementation of the method introduced in this paper is available at https://doi.org/10.5281/zenodo.6983162.

Author contributions

FT and AD designed the method. FT performed the simulations of the method for the mean-field models and analyzed the data. NT-C performed the simulations of the method for LFP in spiking networks. DD and MC designed and performed the simulations for the two-compartment spiking AdEx model. FT, DD, and AD wrote the manuscript. AD supervised the project. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the CNRS and the European Union (Human Brain Project H2020-785907 and H2020-945539).

Acknowledgments

Large scale simulation of the human brain in the Virtual Brain platform are based on previous work of Jennifer S. Goldman at the Paris-Saclay Institute of Neuroscience (NeuroPSI). We thank Alexis Garcia for valuable discussion.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2022.968278/full#supplementary-material

References

Bédard, C., and Destexhe, A. (2012). “Modeling local field potentials and their interaction with the extracellular medium,” in Handbook of Neural Activity Measurement (Cambridge), 136–191.

Google Scholar

Bédard, C., Kröger, H., and Destexhe, A. (2004). Modeling extracellular field potentials and the frequency-filtering properties of extracellular space. Biophys. J. 86, 1829–1842. doi: 10.1016/S0006-3495(04)74250-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Brette, R., and Gerstner, W. (2005). Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J. Neurophysiol. 94, 3637–3642. doi: 10.1152/jn.00686.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., Anastassiou, C. A., and Koch, C. (2012). The origin of extracellular fields and currents–EEG, ECOG, LFP and spikes. Nat. Rev. Neurosci. 13, 407–420. doi: 10.1038/nrn3241

PubMed Abstract | CrossRef Full Text | Google Scholar

De Munck, J., Wolters, C. H., and Clerc, M. (2012). Eeg and meg: forward modeling. Handbook Neural Activity Meas. 19, 192–248. doi: 10.1017/CBO9780511979958.006

CrossRef Full Text | Google Scholar

De Pasquale, F., Della Penna, S., Snyder, A. Z., Lewis, C., Mantini, D., Marzetti, L., et al. (2010). Temporal dynamics of spontaneous meg activity in brain networks. Proc. Natl. Acad. Sci. U.S.A. 107, 6040–6045. doi: 10.1073/pnas.0913863107

PubMed Abstract | CrossRef Full Text | Google Scholar

Depannemaecker, D., Destexhe, A., Jirsa, V., and Bernard, C. (2021). Modeling seizures: from single neurons to networks. Seizure 90, 4–8. doi: 10.1016/j.seizure.2021.06.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., Contreras, D., and Steriade, M. (1999). Spatiotemporal analysis of local field potentials and unit discharges in cat cerebral cortex during natural wake and sleep states. J. Neurosci. 19, 4595–4608. doi: 10.1523/JNEUROSCI.19-11-04595.1999

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Volo, M., Romagnoni, A., Capone, C., and Destexhe, A. (2019). Biologically realistic mean-field models of conductance-based networks of spiking neurons with adaptation. Neural Comput. 31, 653–680. doi: 10.1162/neco_a_01173

PubMed Abstract | CrossRef Full Text | Google Scholar

Garcia-Rodriguez, A., and Destexhe, A. (2021). A formalism for the long-distance magnetic field generated by populations of neurons. arXiv preprint arXiv:2112.02024. doi: 10.48550/arXiv.2112.02024

CrossRef Full Text | Google Scholar

Geselowitz, D. (1970). On the magnetic field generated outside an inhomogeneous volume conductor by internal current sources. IEEE Trans. Magn. 6, 346–347. doi: 10.1109/TMAG.1970.1066765

CrossRef Full Text | Google Scholar

Goldman, J. S., Kusch, L., Yalcinkaya, B. H., Depannemaecker, D., Nghiem, T.-A. E., Jirsa, V., et al. (2021). A comprehensive neural simulation of slow-wave sleep and highly responsive wakefulness dynamics. bioRxiv. doi: 10.1101/2021.08.31.458365

CrossRef Full Text | Google Scholar

Goldman, J. S., Kusch, L., Yalcinkaya, B. H., Depannemaecker, D., Nghiem, T. -A. E., Jirsa, V., et al. (2020). Brain-scale emergence of slow-wave synchrony and highly responsive asynchronous states based on biologically realistic population models simulated in the virtual brain. BioRxiv. doi: 10.1101/2020.12.28.424574

CrossRef Full Text | Google Scholar

Hagen, E., Dahmen, D., Stavrinou, M. L., Lindén, H., Tetzlaff, T., Van Albada, S. J., et al. (2016). Hybrid scheme for modeling local field potentials from point-neuron networks. Cereb. Cortex 26, 4461–4496. doi: 10.1093/cercor/bhw237

PubMed Abstract | CrossRef Full Text | Google Scholar

Hagen, E., Næss, S., Ness, T. V., and Einevoll, G. T. (2018). Multimodal modeling of neural network activity: Computing LFP, ECOG, EEG, and MEG signals with lfpy 2.0. Front Neuroinform. 12, 92. doi: 10.3389/fninf.2018.00092

PubMed Abstract | CrossRef Full Text | Google Scholar

Hämäläinen, M., Hari, R., Ilmoniemi, R. J., Knuutila, J., and Lounasmaa, O. V. (1993). Magnetoencephalography–theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Mod. Phys. 65, 413. doi: 10.1103/RevModPhys.65.413

CrossRef Full Text | Google Scholar

Hansen, P., Kringelbach, M., and Salmelin, R. (2010). MEG: An Introduction To Methods. Oxford: Oxford University Press.

Google Scholar

Ilmoniemi, R. J., and Sarvas, J. (2019). Brain Signals: Physics and mathematics of MEG and EEG. Cambridge: MIT Press.

Google Scholar

Kuhn, A., Aertsen, A., and Rotter, S. (2004). Neuronal integration of synaptic input in the fluctuation-driven regime. J. Neurosci. 24, 2345–2356. doi: 10.1523/JNEUROSCI.3349-03.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindén, H., Hagen, E., Leski, S., Norheim, E. S., Pettersen, K. H., and Einevoll, G. T. (2014). Lfpy: a tool for biophysical simulation of extracellular potentials generated by detailed model neurons. Front. Neuroinform. 7, 41. doi: 10.3389/fninf.2013.00041

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes da Silva, F. H. (2013). EEG and MEG: relevance to neuroscience. Neuron 80, 1112–1128. doi: 10.1016/j.neuron.2013.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandija, S., Petrov, P. I., Vink, J. J., Neggers, S. F., and van den Berg, C. A. (2021). Brain tissue conductivity measurements with mr-electrical properties tomography: an in vivo study. Brain Topogr. 34, 56–63. doi: 10.1007/s10548-020-00813-1

PubMed Abstract | CrossRef Full Text | Google Scholar

McCann, H., Pisano, G., and Beltrachini, L. (2019). Variation in reported human head tissue electrical conductivity values. Brain Topogr. 32, 825–858. doi: 10.1007/s10548-019-00710-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Megıas, M., Emri, Z., Freund, T., and Gulyas, A. (2001). Total number and distribution of inhibitory and excitatory synapses on hippocampal ca1 pyramidal cells. Neuroscience 102, 527–540. doi: 10.1016/S0306-4522(00)00496-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Melozzi, F., Woodman, M. M., Jirsa, V. K., and Bernard, C. (2017). The virtual mouse brain: a computational neuroinformatics platform to study whole mouse brain dynamics. eNeuro 4, ENEURO.0111-17.2017. doi: 10.1523/ENEURO.0111-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Mosher, J. C., Leahy, R. M., and Lewis, P. S. (1999). EEG and MEG: forward solutions for inverse methods. IEEE Trans. Biomed. Eng. 46, 245–259. doi: 10.1109/10.748978

PubMed Abstract | CrossRef Full Text | Google Scholar

Niedermeyer, E., and Lopes da Silva, F. H. (2005). Electroencephalography: Basic Principles, Clinical Applications, and Related Fields. Baltimore, MD: Lippincott Williams and Wilkins.

PubMed Abstract | Google Scholar

Nunez, P. L.Srinivasan, R., et al. (2006). Electric Fields of the Brain: The Neurophysics of EEG. Oxford: Oxford University Press, USA.

PubMed Abstract | Google Scholar

Sanz Leon, P., Knock, S. A., Woodman, M. M., Domide, L., Mersmann, J., McIntosh, A. R., et al. (2013). The virtual brain: a simulator of primate brain network dynamics. Front. Neuroinform. 7, 10. doi: 10.3389/fninf.2013.00010

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanz-Leon, P., Knock, S. A., Spiegler, A., and Jirsa, V. K. (2015). Mathematical framework for large-scale brain network modeling in the virtual brain. Neuroimage 111, 385–430. doi: 10.1016/j.neuroimage.2015.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Schirner, M., Rothmeier, S., Jirsa, V. K., McIntosh, A. R., and Ritter, P. (2015). An automated pipeline for constructing personalized virtual brains from multimodal neuroimaging data. Neuroimage 117, 343–357. doi: 10.1016/j.neuroimage.2015.03.055

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, N. R., Manshanden, I., and da Silva, F. H. L. (2000). A MEG study of sleep. Brain Res. 860, 64–76. doi: 10.1016/S0006-8993(00)01974-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Skaar, J.-E. W., Stasik, A. J., Hagen, E., Ness, T. V., and Einevoll, G. T. (2020). Estimation of neural network model parameters from local field potentials (lfps). PLoS Comput. Biol. 16, e1007725. doi: 10.1371/journal.pcbi.1007725

PubMed Abstract | CrossRef Full Text | Google Scholar

Stokes, M. G., Wolff, M. J., and Spaak, E. (2015). Decoding rich spatial information with high temporal resolution. Trends Cogn. Sci. 19, 636–638. doi: 10.1016/j.tics.2015.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Teleńczuk, B., Dehghani, N., Le Van Quyen, M., Cash, S. S., Halgren, E., Hatsopoulos, N. G., et al. (2017). Local field potentials primarily reflect inhibitory neuron activity in human and monkey cortex. Sci. Rep. 7, 1–10. doi: 10.1038/srep40211

PubMed Abstract | CrossRef Full Text | Google Scholar

Telenczuk, B., Telenczuk, M., and Destexhe, A. (2020). A kernel-based method to calculate local field potentials from networks of spiking neurons. J. Neurosci. Methods 344, 108871. doi: 10.1016/j.jneumeth.2020.108871

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: magnetoencephalography (MEG), local field potential (LFP), mean field (MF), whole brain analysis, forward modeling

Citation: Tesler F, Tort-Colet N, Depannemaecker D, Carlu M and Destexhe A (2022) Mean-field based framework for forward modeling of LFP and MEG signals. Front. Comput. Neurosci. 16:968278. doi: 10.3389/fncom.2022.968278

Received: 13 June 2022; Accepted: 17 August 2022;
Published: 13 October 2022.

Edited by:

Maurizio Mattia, National Institute of Health (ISS), Italy

Reviewed by:

Torbjørn Vefferstad Ness, Norwegian University of Life Sciences, Norway
Ali Khadem, K. N. Toosi University of Technology, Iran

Copyright © 2022 Tesler, Tort-Colet, Depannemaecker, Carlu and Destexhe. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Federico Tesler, ftesler@gmail.com

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.