Numerical investigation of a neural field model including dendritic processing

Daniele Avitabile; Stephen Coombes; Pedro M. Lima

doi:10.3934/jcd.2020011

Article Contents

2020, Volume 7, Issue 2: 271-290. Doi: 10.3934/jcd.2020011

This issue Previous Article Continuous approximation of $ M_t/M_t/ 1 $ distributions with application to production Next Article An energy dissipative semi-discrete finite-difference method on staggered meshes for the 3D compressible isothermal Navier–Stokes–Cahn–Hilliard equations

Numerical investigation of a neural field model including dendritic processing

1.
Department of Mathematics, Vrije Universiteit Amsterdam, De Boelelaan 1081a 1081 HV Amsterdam, The Netherlands
2.
Mathneuro Team, Inria Sophia Antipolis, 2004 Rue des Lucioles, 06902 Sophia Antipolis, Cedex
3.
Centre for Mathematical Medicine and Biology, School of Mathematical Sciences, University of Nottingham, NG7 2RD, UK
4.
CEMAT, Instituto Superior Tecnico, University of Lisbon, Portugal

Received: September 30, 2019

Published: July 2020

Abstract Full Text(HTML) Figure(4) / Table(4) Related Papers Cited by

Abstract

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.

Keywords:

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

Citation:

FullText(HTML)

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

Download: Full-size image PowerPoint slide

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 $

Download: Full-size image PowerPoint slide

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)

Download: Full-size image PowerPoint slide

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 $

Download: Full-size image PowerPoint slide

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

| Show Table

DownLoad: CSV

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

| Show Table

DownLoad: CSV

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

| Show Table

DownLoad: CSV

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 $

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

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