\`x^2+y_1+z_12^34\`
Advanced Search
Article Contents
Article Contents

Numerical investigation of a neural field model including dendritic processing

Abstract Full Text(HTML) Figure(4) / Table(4) Related Papers Cited by
  • 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:

    \begin{equation} \\ \end{equation}
  • 加载中
  • 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
     | 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
  • [1] U. M. AscherS. 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. BojakT. F. OostendorpA. 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. CrookG. B. ErmentroutM. 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. HeitmannM. 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. VisserR. NicksO. 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)

SHARE

Article Metrics

HTML views(478) PDF downloads(262) Cited by(0)

Access History

Other Articles By Authors

Catalog

    /

    DownLoad:  Full-Size Img  PowerPoint
    Return
    Return