Generalized cable formalism to calculate the magnetic field of single neurons and neuronal populations

Neurons generate magnetic fields which can be recorded with macroscopic techniques such as magnetoencephalography. The theory that accounts for the genesis of neuronal magnetic fields involves dendritic cable structures in homogeneous resistive extracellular media. Here we generalize this model by considering dendritic cables in extracellular media with arbitrarily complex electric properties. This method is based on a multiscale mean-field theory where the neuron is considered in interaction with a “mean” extracellular medium (characterized by a specific impedance). We first show that, as expected, the generalized cable equation and the standard cable generate magnetic fields that mostly depend on the axial current in the cable, with a moderate contribution of extracellular currents. Less expected, we also show that the nature of the extracellular and intracellular media influence the axial current, and thus also influence neuronal magnetic fields. We illustrate these properties by numerical simulations and suggest experiments to test these findings.

Figure 1

Scheme to calculate the magnetic induction produced by a dendritic branch. (a) To evaluate the contribution of the dendritic segment, we divide space into three regions: L, P, R. We first evaluate ${\mathbf{B}}^{\theta }$ in the principal region P, which corresponds to the space between Regions L and R. Next, we evaluate ${\mathbf{B}}^{\theta }$ in the boundary regions $L$ and $R$. Note that the knowledge of ${\mathbf{B}}^{\theta }$ in Region $P$ is necessary to evaluate ${\mathbf{B}}^{\theta }$ in Regions $L$ and $R$ because one must know ${\mathbf{B}}^{\theta }$ on the two planes $z=0$ and $z={\sum }_{i=1}^{N_p}{l}_i=l$, in order to calculate its explicit value in Regions $L$ and $R$ using Eq. (5). (b) Evaluation of ${\mathbf{B}}^{\theta }$ for a segment of variable diameter. In this case, the same procedure is followed, except that Region P is divided into $N_p$ compartments, each described by a continuous cylinder, $P={\bigcup }_{i=1}^{N_p}{p}_i$. Note that the continuity conditions on the axial current and the transmembrane voltage allow one to define boundary conditions for ${\mathbf{B}}^{\theta }$ over the surfaces of the compartments $p_i$. The figure shows an example with $N_p=3$.

Figure 2

Coordinate scheme for a cable segment of constant diameter. The scheme shows the cable with the cylindric coordinate system used in the paper, as well as the surfaces $A$ and $C$, which are the sections that cuts the cable perpendicular to its membrane (delimited by surface $B)$. $\mathcal{D}$ is the interior volume of the segment, as delimited by these surfaces, and $\partial \mathcal{D}$ is the reunion of the two surfaces $A$ and $B$.

Figure 3

Example with two neurons. In order to calculate the value of the magnetic induction $\vec{\mathbf{B}}$ generated by many neurons, one has to sum the values of ${\vec{\mathbf{B}}}_i$ produced by each branch. Thus, it is sufficient to know the axial current $i_i^g$ at each branch to calculate $\vec{\mathbf{B}}$.

Figure 4

Illustration of the current fields around the soma of a ball-and-stick model. The current fields are shown (arrows) around the soma when the generalized membrane current is perpendicular to the soma membrane (red arrows). The isopotential surfaces are shown in blue and correspond to the soma membrane. If the soma has a different “diameter” but coincides with the isopotential surface, then the geometry of these current lines and isopotential surfaces remains invariant. However, the value of the electric potential is different on each equipotential surface.

Figure 5

Equivalent scheme to calculate the current flowing from distal to proximal at the position of the synapse, when the synaptic current is known.

Figure 6

Synaptic current sources used in the simulations. (a) Example of excitatory (blue, top curve) and inhibitory (black, bottom curve) current sources used in simulations. These examples consists of 1000 random synaptic events per second. (b) and (c) Modulus and phase, respectively, of the complex Fourier transform of these processes. Note that the inhibitory current in not represented in (b) because its modulus is identical to that of the excitatory current. The red dashed line in (b) corresponds to a Lorentzian $\left(\frac{A}{1+i\omega {\tau }_m}\right)$ with ${\tau }_m=5\mathrm{ms}$ and $\left|A\right|=1\mathrm{nA})$.

Figure 7

Magnetic induction for the resistive model. ${\mathbf{B}}^{\theta }$ is shown here at the surface of the dendrite, as a function of position (distance to soma) for different frequencies between 1 Hz and $5000\mathrm{Hz}$. The blue dashed lines correspond to ${\mathbf{B}}^{\theta }$ generated when only excitatory synapses were present, and the black curves correspond to both synapses present. ${\mathbf{B}}^{\theta }$ is always decreasing with frequency and is larger and approximately constant between the two locations of the synaptic currents.

Figure 8

Frequency profile of the magnetic field for the resistive model. ${\mathbf{B}}^{\theta }$ at the surface of the dendrite is represented as a function of frequency at different positions (both excitatory and inhibitory synapses were present). The red curves correspond to different positions between the inhibitory synapses and the soma, the blue curves are taken at different positions between the excitatory synapses and the end of the dendrite, and the black curves represent positions in between the two synapse sites. Note that the modulus of ${\mathbf{B}}^{\theta }$ does not depend on position.

Figure 9

Magnetic induction ${\mathbf{B}}^{\theta }$ on the surface of the dendrite for a neuron embedded in diffusive media. ${\mathbf{B}}^{\theta }$ is represented for different frequencies. The blue dashed curves correspond to ${\mathbf{B}}^{\theta }$ produced at the surface of the dendrite with only excitatory synapses, and black curves correspond to excitatory and inhibitory synapses present. We see that ${\mathbf{B}}^{\theta }$ is a decreasing function of frequency and is higher towards inhibitory synapses, and low outside of this region.

Figure 10

Magnetic induction ${\mathbf{B}}^{\theta }$ on the surface of the dendrite, as a function of frequency, for a neuron within diffusive media. ${\mathbf{B}}^{\theta }$ is represented for different positions on the dendrite, with both excitatory and inhibitory synapses present. One can see three distinct regions: proximal region between the soma and the location of inhibitory synapses (red curves), region between the two synaptic sites (black curves), and the distal region between the location of excitatory synapses and the end of the dendrite (blue curves). Note that the modulus of ${\mathbf{B}}^{\theta }$ depends very weakly on dendritic position when we are in between the two synaptic sites.

Figure 11

Distance dependence of the magnetic induction for a ball-and-stick model with resistive media. The boundary conditions are represented in Figs. 10 and 11. (a) Attenuation law for the modulus of ${\mathbf{B}}^{\theta }$ relative to $r$ (direction perpendicular to the axis of the stick). For $r<100\mu \mathrm{m}=l/6$, the attenuation is varying as $1/r$ with a proportionality constant that depends on frequency. For $r>200\mu \mathrm{m}=l/3$, the attenuation varies as $1/r^2$ and is roughly independent of frequency. (b) Attenuation law relative to $z$ (direction parallel to the axis of the stick). The attenuation does not depend on frequency for positions outside the regions between the synapses. In all cases, the phase varied very little and was not represented.

Figure 12

Distance dependence of the magnetic induction for a ball-and-stick model with diffusive media. The same arrangement as in Fig. 11, but with boundary conditions as represented in Figs. 9 and 10. (a) Attenuation law for the modulus of ${\mathbf{B}}^{\theta }$ relative to $r$. As for the resistive model, the attenuation varies as $1/r$ for $r<100\mu \mathrm{m}=l/6$, and as $1/r^2$ for $r>200\mu \mathrm{m}=l/3$. (b) Attenuation law relative to $z$. Contrary to the resistive model, the attenuation depends on frequency for all positions.

Figure 13

Scheme to calculate ${\mathbf{B}}^{\theta }$. (a) Extension of the compartment in Regions L and R when the principal region (P) consists of a single cylinder compartment of radius $a$. (b) Calculation scheme. In a first step, we calculate in Fourier space the field ${\mathbf{B}}^{\theta }$ by assuming that the boundary conditions on the cylinder are such that ${\mathbf{B}}^{\theta }(a,z<0,\omega )=0$ and ${\mathbf{B}}^{\theta }(a,z>l,\omega )=0$ in Regions L and R. The solution of ${\nabla }^2\vec{\mathbf{B}}=0$ obtained in such conditions is called the first-order solution of ${\mathbf{B}}^{\theta }$. In a second step, we improve the boundary conditions by applying the Hankel transform of order 1 relative to $r$, and the continuity principle of the solution at the borders L-P and P-R. We reevaluate the solution of ${\nabla }^2\vec{\mathbf{B}}=0$ with these new boundary conditions to obtain a second-order solution of ${\mathbf{B}}^{\theta }$. The same iteration is continued.

Figure 14

Bessel functions for a continuous cylinder compartment. The Bessel functions are indicated as a function of the distance $r$ perpendicular to the axis of the cylinder; the cylinder had a length $l$ of $300\mu \mathrm{m}$ and a radius $a=2\mu \mathrm{m}$. We have $k^{\prime }=2\times {10}^4{\mathrm{m}}^{-1}$. (a) $K_1$ as a function of distance. The red curve shows the function $K_1\left(kr\right)$ with $k=k^{\prime }$. The black dashed straight line represents the function $1/kr$, and the blue dashed curve represents the asymptotic behavior of $K_1$ for $r\to \infty$. We have $K_1\left(kr\right)\stackrel{\infty }{}\sqrt{\frac{\pi }{2kr}}{e}^{-kr}$ . At short distances (smaller than $2/k)$, the function $K_1$ decays linearly with distance, but for large distances $\left(r>\pi /k\right)$, it converges more rapidly than an exponential decay with distance. (b) The function $I_1$ (modified Bessel function of first kind of order 1) is well approximated by a straight line $\left[I_1\left(kr\right)=kr/2\right]$ for $k=k^{\prime }$ when $r<a$. (c) Bessel function of first kind of order 1 when $r>a$. The blue and black curves correspond, respectively, to $k=k^{\prime }$ and $k=5k^{\prime }$. Note that we have $J_1\left(kr\right)\stackrel{\infty }{}\sqrt{\frac{2}{\pi kr}}cos(kr-\frac{3\pi }{4})$. We see that when $k$ is large enough, the function $J_1\left(kr\right)$ can capture small spatial variations. (d) The Bessel function of first kind of order 1 is equivalent to a straight line $\left[J_1\left(kr\right)=kr/2\right]$ when $r<a$. The blue dashed curve corresponds to approximating $J_1\left(kr\right)$ by a linear law for $k=k^{\prime }$, while the black dashed curve is the linear approximation for $k=5k^{\prime }$. The red curves correspond to $J_1\left(kr\right)$ for $k=k^{\prime }$ and $k=5k^{\prime }$ [33, 34].

Figure 15

Example of application of the Hankel transform of order 1. (a) Approximation using the Hankel transform of order 1 of the function $f\left(r\right)=H(r-a)1/r$ with $a=1\mu \mathrm{m},1<k_r<5\times {10}^6$, and $\Delta k_r={10}^3$. The values smaller than ${10}^3$ are not significant because of the value of $\Delta k_r$ is larger than 1000. We can see that the approximation using the Hankel transform is valid for distances up to $1\mathrm{mm}$. (b) Inverse transform applied to this approximation (in blue), and comparison with the original function (in red), between $1\mu \mathrm{m}$ and $1\mathrm{mm}$. The parameter $k_r$ of the Hankel transform plays a similar role as the wave number $\left(\frac{2\pi }{\lambda }\right)$ in spatial Fourier transform. The larger $k_r$, the more sensitive to fine spatial details.

Figure 16

Organigram of the iterative method to calculate ${\mathbf{B}}^{\theta }$.