Article Contents
Article Contents

# Numerical investigation of a neural field model including dendritic processing

• We consider a simple neural field model in which the state variable is dendritic voltage, and in which somas form a continuous one-dimensional layer. This neural field model with dendritic processing is formulated as an integro-differential equation. We introduce a computational method for approximating solutions to this nonlocal model, and use it to perform numerical simulations for neuro-biologically realistic choices of anatomical connectivity and nonlinear firing rate function. For the time discretisation we adopt an Implicit-Explicit (IMEX) scheme; the space discretisation is based on a finite-difference scheme to approximate the diffusion term and uses the trapezoidal rule to approximate integrals describing the nonlocal interactions in the model. We prove that the scheme is of first-order in time and second order in space, and can be efficiently implemented if the factorisation of a small, banded matrix is precomputed. By way of validation we compare the outputs of a numerical realisation to theoretical predictions for the onset of a Turing pattern, and to the speed and shape of a travelling front for a specific choice of Heaviside firing rate. We find that theory and numerical simulations are in excellent agreement.

Mathematics Subject Classification: Primary: 58F15, 58F17; Secondary: 53C35.

 Citation:

• Figure 1.  Schematic of the neural field model. (a) Dendrites are represented as unbranched fibres (red), orthogonal to a continuum of somas (somatic layer, in grey). (b) Model (1) is for a 1D somatic layer, with coordinate $x \in \mathbb{R}$, and fiber coordinate $\xi \in \mathbb{R}$. Input currents are generated in a small neighbourhood of the somatic layer, at $\xi = 0$ and are delivered to a contact point, in a small neighbourhood of $\xi = \xi_0$. The strength of interaction depends on the distance between sources and contact points, measured along the somatic layer, hence the inputs that are generated at $A$ and transmitted to $B$, $C$, and $D$ depend on $|x_B - x_A|$, $|x_C-x_A|$, and $|x_D - x_A|$, respectively (see (2))

Figure 2.  Coherent structure observed in time simulation of (3), (22), (23). (a): Pseudocolor plot of $V(x,\xi,t)$ at several time points, showing two counter-propagating waves. (b): Solution at $\xi = 0$, showing the wave profile. Parameters: $\xi_0 = 1$, $\varepsilon = 0.005$, $\nu = 0.4$, $\gamma = 1$, $\beta = 1000$, $\theta = 0.01$, $\kappa = 3$, $L_x = 24 \pi$, $L_\xi = 3$, $n_x = 2^{10}$, $n_\xi = 2^{12}$, $\tau = 0.05$

Figure 3.  (a) Travelling wave speed versus firing rate threshold, computed analytically via (24), and numerically via the time-stepper. (b)-(d) Convergence of the computed speed to the analytical speed at $\theta = 0.01$, as a function of the the kernel support parameter $\varepsilon$, the steepness of the sigmoid $\beta$, and the time-stepping parameter $\tau$, respectively. Parameters as in (2)

Figure 4.  Numerical simulation of Turing bifurcation for the model with kernel and firing rate function given in (28). (a): Plots of $\hat w(p) - w_*$ for $\beta = 28$ and $\beta = 30$, from which we deduce that a Turing bifurcation occurs for an intermediate value of $\beta$. We expect perturbations of the trivial state to decay exponentially if $\beta = 28$ and to increase exponentially if $\beta = 30$, as confirmed in panels (b)-(d). (b): Maximum absolute value of the voltage at $\xi = 0$, as function of time, when the trivial steady state is perturbed. (c), (d): Pseudocolor plot of $V$ when $\beta = 28$ and $\beta = 30$, respectively. Parameters: $\xi_0 = 1$, $\varepsilon = 0.005$, $\nu = 6$, $c = 1$, $a_1 = 1$, $b_1 = 1$, $a_2 = 1/4$, $b_2 = 1/2$, $n_x = 2^9$, $L_x = 10\pi$, $n_\xi = 2^{11}$, $L_\xi = 2.5 \pi$, $\tau = 0.01$

 Algorithm 1: IMEX time stepper in matrix form (8), nonlinear term computed with pseudospectral evaluation (21) Input: Initial condition $V^0 \in \mathbb{R}^{n_\xi \times n_x}$, time step $\tau$, number of steps $n_t$. Output: An approximate solution $(V^n)_{n=1}^{n_t} \subset \mathbb{R}^{n_\xi \times n_x}$ 1 begin 2      Compute grid vectors $\xi \in \mathbb{R}^{n_\xi \times 1}$, $x \in \mathbb{R}^{1 \times n_x}$. 3       Compute synaptic vectors $w, \hat w = \mathcal{F}_{n_x}[w] \in \mathbb{R}^{1 \times n_x}$. 4       Compute synaptic vectors $\alpha, \alpha' \in \mathbb{R}^{n_\xi \times 1}$. 5       Compute quadrature weights $\sigma \in \mathbb{R}^{n_\xi \times 1}$. 6       Compute sparse $LU$-factorisation of $A$, $LU = A \in \mathbb{R}^{n_\xi \times n_\xi}.$ 7       for $n = 1,\ldots,n_t$ do 8              Set $V = V^{n-1} \in \mathbb{R}^{n_\xi \times n_x}$. 9               Compute the external input at time $t_{n-1}$ and store it in $G \in \mathbb{R}^{n_\xi \times n_x}$. 10               Set $z = \mathcal{F}_{n_x}\big[(\alpha' \odot \sigma)^T S(V)\big] \in \mathbb{R}^{1 \times n_x}$. 11                Set $N = h_x \alpha \mathcal{F}^{-1}_{n_x} [ \hat w \odot z] \in \mathbb{R}^{n_\xi \times n_x}$. 12                Solve for $V^{n}$ the linear problem $(LU)V^n = V + \tau(N+G)$. 13      end 14 end
 Algorithm 2: IMEX time stepper in vector form (7), nonlinear term evaluated with quadrature formula (12). Input: Initial condition $U^0 \in \mathbb{R}^{n_\xi n_x}$, time step $\tau$, number of steps $n_t$. Output: An approximate solution $(U^n)_{n=1}^{n_t} \subset \mathbb{R}^{n_\xi n_x}$ 1 begin 2      Compute grid vectors $\xi \in \mathbb{R}^{n_\xi}$, $x \in \mathbb{R}^{n_x}$. 3       Compute synaptic vector $w \in \mathbb{R}^{1 \times n_x}$. 4       Compute synaptic vectors $\alpha, \alpha' \in \mathbb{R}^{n_\xi \times 1}$. 5       Compute quadrature weights $\rho \in \mathbb{R}^{n_x}$, $\sigma \in \mathbb{R}^{n_\xi}$. 6      Compute sparse $LU$-factorisation of A, $LU = \big( (1+\tau \gamma) I_{n_x n_\xi} - \tau \nu I_{n_x} \otimes D_{\xi \xi} \big) \in \mathbb{R}^{n_\xi n_x \times n_\xi n_x }.$ 7      for $n = 1,\ldots,n_t$ do 8              Set $Z = U^{n-1} \in \mathbb{R}^{n_\xi n_x}$. 9              Compute the external input at time $t_{n-1}$ and store it in $G \in \mathbb{R}^{n_\xi n_x}$. 10              Compute the nonlinear term $N$ using (12). 11              Solve for $U^{n}$ the linear problem $(LU)U^n = Z + \tau(N+G)$. 12      end 13 end

Table 1.  Flop count for the initialisation step (lines 2-6) and for one time step (lines 8-12) in Algorithms 1, 2

 Algorithm 1 Algorithm 2 Lines Flops Lines Flops 2 $n_\xi + n_x$ 2 $n_\xi + n_x$ 3 $2n_x$ 3 $n_x$ 4 $2n_\xi$ 4 $2n_\xi$ 5 $n_\xi$ 5 $n_\xi + n_x$ 6 $2n_\xi-1$ 6 $2n_\xi n_x -1$ 2-6 $O(n_\xi) + O(n_x)$ 2-6 $O(n_\xi n_x)$ 8 $n_\xi n_x$ 8 $n_\xi n_x$ 9 $n_\xi n_x$ 9 $n_\xi n_x$ 10 $3 n_\xi n_x + O(n_x \log n_x) + n_\xi - n_x$ 10 $2n^2_\xi n_x^2 - n^2_\xi n_x$ 11 $n_\xi n_x + O(n_x \log n_x) + 2n_x$ 11 $5n_\xi n_x -4$ 12 $5 n_\xi n_x - 4 n_x$ $8-12$ O(n_\xi n_x) + O(n_x \log n_x) $8-11$ O(n^2_\xi n^2_x) $Table 2. Space requirements, measured in Floating Point Numbers, for Algorithms 1 and 2. Arrays$ d_1, \ldots, d_3 $, store diagonals of the$ LU $-factorisation in the respective algorithms  Floating Point Numbers Algorithm 1 Algorithm 2$ n_\xi  \xi,\alpha,\alpha',d_1,d_2,d_3,z  \xi, \rho, \alpha, \alpha'  n_x  x,w,\sigma  x, \sigma, w  n_\xi n_x  V,V^n,G,N  U^n,Z,N,G,d_1,d_2,d_3 $Total$ 4n_x n_\xi + 7n_\xi + 3n_x  7 n_x n_\xi +2n_\xi + 2n_x \$
•  [1] U. M. Ascher, S. J. Ruuth and B. T. R. Wetton, Implicit-explicit methods for time-dependent partial differential equations, SIAM J. Numer. Anal., 32 (1995), 797-823.  doi: 10.1137/0732037. [2] D. Avitabile, S. Coombes and P. M. Lima, danieleavitabile/neural-field-with-dendrites: Ancillary codes to "Numerical Investigation of a Neural Field Model Including Dendritic Processing", 2020. doi: 10.5281/zenodo.3731920. [3] I. Bojak, T. F. Oostendorp, A. T. Reid and R. Kötter, Connecting mean field models of neural activity to EEG and fMRI data, Brain Topography, 23 (2010), 139-149.  doi: 10.1007/s10548-010-0140-3. [4] P. C. Bressloff, New mechanism for neural pattern formation, Phys. Rev. Lett., 76 (1996), 4644-4647.  doi: 10.1103/PhysRevLett.76.4644. [5] P. C. Bressloff and S. Coombes, Physics of the extended neuron, International Journal of Modern Physics B, 11 (1997), 2343-2392.  doi: 10.1142/S0217979297001209. [6] S. Coombes, Waves, bumps, and patterns in neural field theories, Biol. Cybernet., 93 (2005), 91-108.  doi: 10.1007/s00422-005-0574-y. [7] S. Coombes, P. Beim Graben, R. Potthast and J. Wright, Neural Fields: Theory and Applications, Springer, 2014. [8] S. Coombes, H. Schmidt and I. Bojak, Interface dynamics in planar neural field models, J. Math. Neurosci., 2 (2012), 9, 27 pp. doi: 10.1186/2190-8567-2-9. [9] S. M. Crook, G. B. Ermentrout, M. C. Vanier and J. M. Bower, The role of axonal delay in the synchronization of networks of coupled cortical oscillators, Journal of Computational Neuroscience, 4 (1997), 161-172. [10] S. Heitmann, M. J. Aburn and M. Breakspear, The Brain dynamics toolbox for matlab, Neurocomputing, 315 (2018), 82-88.  doi: 10.1016/j.neucom.2018.06.026. [11] A. Hutt and N. Rougier, Numerical simulation scheme of one-and two dimensional neural fields involving space-dependent delays, Neural Fields, Springer, Heidelberg, (2014), 175–185. [12] P. M. Lima and E. Buckwar, Numerical solution of the neural field equation in the two-dimensional case, SIAM J. Sci. Comput., 37 (2015), B962–B979. doi: 10.1137/15M1022562. [13] E. J. Nichols and A. Hutt, Neural field simulator: Two-dimensional spatio-temporal dynamics involving finite transmission speed, Front. Neuroinform., 9 (2015), 25. doi: 10.3389/fninf.2015.00025. [14] P. L. Nunez, Neocortical dynamics and human EEG rhythms, Physics Today, 49 (1996), 57. doi: 10.1063/1.2807585. [15] K. H. Pettersen and G. T. Einevoll, Amplitude variability and extracellular low-pass filtering of neuronal spikes, Biophysical Journal, 94 (2008), 784-802.  doi: 10.1529/biophysj.107.111179. [16] J. Rankin, D. Avitabile, J. Baladron, G. Faye and D. J. B. Lloyd, Continuation of localized coherent structures in nonlocal neural field equations, SIAM J. Sci. Comput., 36 (2014), B70–B93. doi: 10.1137/130918721. [17] J. Ross, M. Margetts, I. Bojak, R. Nicks, D. Avitabile and S. Coombes, A brain-wave equation incorporating axo-dendritic connectivity, Physical Review E, 101 (2020), 022411. [18] P. Sanz-Leon, P. A. Robinson, S. A. Knock, P. M. Drysdale, R. G. Abeysuriya, F. K. Fung, C. J. Rennie and X. Zhao, NFTsim: Theory and simulation of multiscale neural field dynamics, PLoS Computational Biology, 14 (2018), e1006387. doi: 10.1371/journal.pcbi.1006387. [19] L. N. Trefethen, Spectral Methods in MATLAB, Software, Environments, and Tools, SIAM, vol. 10, 2000. doi: 10.1137/1.9780898719598. [20] J. M. Varah, A lower bound for the smallest singular value of a matrix, Linear Algebra and Applications, 11 (1975), 3-5.  doi: 10.1016/0024-3795(75)90112-3. [21] S. Visser, R. Nicks, O. Faugeras and S. Coombes, Standing and travelling waves in a spherical brain model: The Nunez model revisited, Physica D, 349 (2017), 27-45.  doi: 10.1016/j.physd.2017.02.017.

Figures(4)

Tables(4)