Neural field simulator: two-dimensional spatio-temporal dynamics involving finite transmission speed

Neural Field models (NFM) play an important role in the understanding of neural population dynamics on a mesoscopic spatial and temporal scale. Their numerical simulation is an essential element in the analysis of their spatio-temporal dynamics. The simulation tool described in this work considers scalar spatially homogeneous neural fields taking into account a finite axonal transmission speed and synaptic temporal derivatives of first and second order. A text-based interface offers complete control of field parameters and several approaches are used to accelerate simulations. A graphical output utilizes video hardware acceleration to display running output with reduced computational hindrance compared to simulators that are exclusively software-based. Diverse applications of the tool demonstrate breather oscillations, static and dynamic Turing patterns and activity spreading with finite propagation speed. The simulator is open source to allow tailoring of code and this is presented with an extension use case.


INTRODUCTION
The understanding of spatio-temporal electric activity in neural tissue is essential in the study of neurobiological phenomena. To achieve this, mesoscopic-scale models such as neural mass and neural fields (NFM) which describe the dynamics of a large population of neurons reflecting coarse-grained properties of single neurons (Wilson and Cowan, 1973;Deco et al., 2008;Bressloff, 2012;Hutt and Buhry, 2014) play an important role. NFMs serve as a good description of the dynamic source of Local Field Potentials and encephalographic data (Nunez, 1974(Nunez, , 2000Wright and Kydd, 1992;Wright and Liley, 1994;Jirsa et al., 2002;Nunez and Srinivasan, 2006;Coombes et al., 2014). They allow to consider diverse single neuron features that may tune neural population dynamics, such as somatic (Molaee-Ardekani et al., 2007) and synaptic adaptation (Coombes and Owen, 2005;Kilpatrick and Bressloff, 2010), extrasynaptic receptor dynamics (Hashemi et al., 2014;Hutt and Buhry, 2014) or finite axonal transmission speed (Jirsa and Haken, 1996;Pinto and Ermentrout, 2001;Hutt et al., 2003Hutt et al., , 2008Coombes, 2005;Faye and Faugeras, 2010;Faugeras, 2011, 2013). All these applications make NFMs valuable in order to understand spatio-temporal dynamics of neural population activity.
Mathematical analysis and the numerical integration of NFMs are complementary. The recent years have shown strong attention of research on the mathematical properties of NFMs, whereas the numerical simulation of NFM solutions has been less considered in research. Since NFMs generalize partial differential equations Hutt, 2007) while involving finite transmission delay interactions, they allow to study a large class of pattern forming systems, cf. Hutt (2007). In recent years, several software tools have been developed to simulate neural network dynamics. Examples for simulators for networks of spiking neurons are BRIAN (Stimberg et al., 2014) and Neuron (Carnevale and Hines, 2006). The Virtual Brain (Sanz Leon et al., 2013) allows to simulate networks of neural mass models to reproduce global brain activity. The simulation platform DANA (Rougier and Fix, 2012) simulates a hierarchy of coupled Dynamic Neural Fields which are decentralized, i.e., are updated numerically in time asynchronously (Rougier and Hutt, 2011). These latter software tools are powerful, general and highly adaptive to the framework of their neural network types. However, they do not provide the effective computation for the specific NFM given in Equation (1) which is a stochastic delayed integral-differential equation in two spatial dimensions. The tool presented here fills a gap in the landscape of neural simulator tools which are typically very general and adaptive and, hence, not efficient for NFM. A simulation tool for NFM allows to explore rapidly and in a user-friendly way the solution space of Equation (1) in order to reproduce numerically experimental spatio-temporal dynamics, e.g., to understand neuroimaging data Pinotsis and Friston, 2014), retrieve neural sources and lateral connections (Pinotsis et al., 2013), and understand power spectra of electroencephalographic activity (Pinotsis et al., 2012). In addition, the tool promises to allow detection of new numerical solutions, cf. Section 3. The numerical analysis is non-trivial and challenging if NFMs become more complex, e.g., by involving complex dynamical features rendering the model high-dimensional or by considering delayed interactions. The present work considers a twodimensional spatial embedding of neural populations similar to several previous studies (Laing, 2005;Owen et al., 2007) while taking into account finite axonal transmission speed Rougier, 2010, 2014). By virtue of its modularity, the simulator allows subsequent extensions with additional features, such as extra-synaptic receptor effects or several interacting populations.
The combination of finite axonal transmission speed and two-dimensional spatial embedding is challenging from a numerical simulation perspective due to the missing convolution structure Rougier, 2010, 2014) leading to long simulation durations. To overcome this problem, a numerical technique has been developed in recent years Rougier, 2010, 2014). Since future research in neural fields will investigate spatio-temporal dynamics involving finite axonal transmission speed, we have developed an open-source simulation toolbox that allows to gain spatio-temporal solutions of NFM models in two spatial dimensions, visualize them and save them, if necessary, as movies. We hope that the tool will provide an essential tool for the computational neuroscience community to advance the research field and the insight into the brain.
The simulator in this work obeys integral-differential equations of the type with a two-dimensional square spatial domain and periodic boundary condition. The mean neuron potential V ∈ R at location x ∈ is evolved by the external stimulus I(x, t) ∈ R and the integral of the synaptic connectivity kernel K : R 2 → R and population firing rate S ∈ R which depend on the distance between spatial locations x and y with a finite axonal transmission speed c. Equation (1) represents the core of most NFM in the sense that most NFMs consider extensions of this equation.
Motivation for the work arises from a need for a visualization tool that is useful to the largest number of NFM researchers, allows for the tailoring of code and has fast while visually appealing output. The simulator can operate on all major operating systems. Output of data in three dimensions is provided by PyOpenGL which brings the speed and graphical detail of low-level OpenGL to the agile Python language. It is open source, enabling modification of the simulator in any beneficial way.

MATERIALS AND METHODS
The cross-platform simulator is written in Python (version 2.7) and uses the NumPy library in consideration of its speed being close to the computational rate of the platformdependent C language (Langtangen, 2006). The simulator can be downloaded 1 in a package along with documentation 2 describing its installation, running, features, and examples and the code is registered in ModelDB 3 .
The following Section 2.1 describes the comprehensive access to field parameters, the subsequent section details the techniques applied to accelerate the simulation and Section 2.3 discusses the 3D visualization.

Field Parameters
A textual interface named values.py is provided in the root directory of the simulator code. It allows field values to be changed without knowledge of the inner workings of the simulator. For example, if η in Equation (1) is initialized as a non-zero number, a second order derivative is calculated to solve V. Conversely, the interface eliminates the knowledge requirement of the numerical implementation of the derivatives and all other underlying code implementations. The interface has additional benefits of easily modifying variables in one place without searching through the code. This implementation permits changing parameters easily and sharing code amongst others working with similar simulations by the exchange of a single file.
The most important aspect of a text-based interface from its user experience is its facilitation of novelty by allowing absolute control of all terms of Equation (1). For instance, matrix I can be defined in the interface with as many lines of Python code as necessary given the definition ends with an assignment (i.e., I = ...). Assigning the first parameter in the values.py file, named showData, a value of 3 displays I in the simulator, which can be useful when refining its values. Time-varying spatiotemporal input is available in the interface by uncommenting and modifying the body of a function named updateI in any manner while maintaining that I is returned. Neural field investigations are thereby efficiently implemented with free choice over all the variables accessible through the interface while retaining the full performance.

Accelerated Simulation
The simulator is advantageous in its acceleration of spatial and temporal integration. Multiple approaches are used to increase the simulation speed.

Spatial and Temporal Integration
Equation (1) includes a spatial integral with a homogeneous kernel K. Please recall that homogeneous kernels just depend on the difference vector x − y between two spatial locations x, y including isotropic kernels, i.e., K = K( x − y ), as a specific case. In the absence of the finite transmission delay term, this integral would represent a spatial convolution and would be solvable numerically efficiently by a Fast Fourier Transform (FFT) (Van Loan, 1991). For non-vanishing transmission delay, the convolution structure is less obvious and the FFT is not applicable directly. Nevertheless, it is possible to re-write the spatial integration to utilize a FFT in space (Owen et al., 2007;Rougier, 2010, 2014) as with the maximum delay time τ m and the spatio-temporal kernel function L( We observe that the spatial summation represents an integration over delayed spatial rings, which are convolved spatially with the transfer function S in Equation (1). Introducing a regular rectangular spatial grid for spatial discretization, finite axonal speed c yields rings of width delineated within the field, where n and l are the number of discretized spatial units and the length of the field, respectively, and t is the finite integration time step. The Pythagorean theorem gives the maximum radius of the rings in the field r = n/ √ 2 over which the spatial integration is performed, which is applied to obtain the number of rings in a field as defining the maximum delay to τ m = n rings t. The spatiotemporal kernel L is determined by the spatial kernel K and the axonal speed c Rougier, 2010, 2014). Equation (1) involves distance-dependent delays which represent a specific type of distributed delays (Hutt and Lefebvre, in press). To this end, it is necessary to initialize the field variable V in an initial time interval and the toolbox allows the user to set the initial values arbitrarily. The external input I may be deterministic or stochastic and the user may choose it according to her needs, e.g., implementing spatial correlations in stochastic inputs. To integrate the evolution equation in time the user may choose between different integration methods for delay differential equations Winkler, 2006, 2007). Standard methods discretize time regularly in steps of duration t yielding results (Equations 2, 3). In the case of stochastic input, the toolbox includes numerical implementations of the delayed Euler-Maruyama method (Buckwar et al., 2008) and the stochastic version of the Runge-Kutta method for delayed differential equation (Carletti, 2006). For deterministic inputs, the equivalent deterministic methods are available.
If there is no modification to K and c during the simulation, then L is calculated once only before the start of the simulation while S(·) changes over time. The convolution of L and S is performed using a FFT what greatly increases the speed of the integral convolution compared to conventional integration. This can be understood easily recalling that the two-dimensional FFT needs to sum up n 2 log 2 2 (n) terms leading to summands of the total number of N FFT = n 2 log 2 2 (n)n rings . In contrast, conventional integration sums up terms of number n 4 for each delay time and hence the total number of computation N conv = n 4 n rings , cf. Appendix I. Hence the FFT implementation speeds up the integration by a factor of The axonal speed's implementation is described in detail in Hutt and Rougier (2010). It is important to note that other (conventional) numerical software tools not taking into account the convolution structure have to sum up N conv terms in case of fully connected networks. For instance, this holds true for the simulation tools BRIAN, Neuron and The Virtual Brain (Sanz Leon et al., 2015) which have to memorize the history and advance the stored activity field n rings times. We also note that the method proposed may be implemented in these simulation tools in the future since they also may consider spatially homogeneous neural fields as specific cases. The FFT-based method presented here computes the network interactions faster than these tools by f speedup given in Equation (4). For instance, for a typical number of spatial grid intervals of n = 512 as used in the application showed in Section 3, we obtain the huge speed up factor of f speedup =≈ 3236.

Self-writing Code
The second approach to increase the simulation speed employs self-writing code to reduce the simulator's instruction set. The simulator writes and executes its own code to increase the efficiency of simulations and display only the user defined features. The simulation code is based on interface selections and is self-written by an initialization module at the onset of the program. The interface offers features, such as a second derivative calculation, I and K updates and added noise, that conditionally run during the simulation and are not performed over time if the user selects to view V, I, or K at t = 0. For example, the visual interface offers the choice of viewing the spatial kernel K for its design and visualization. Only the code that initializes and displays K is written to the executing module if this choice is selected. The self-writing code is also favorable when the full simulation is run with calculations executed unconditionally. The result is very efficient code that is changed with every modification to the interface.

Implementation on GPUs
The third approach parallelizes the output calculations on the graphics processing unit (GPU). The displayed matrix is put onto the running system's GPU for hardware acceleration of the visualization. Vertex buffer objects improve visualization throughput by uploading vertices to video device memory where vertex and fragment shaders transform and write neural field data in parallel to the framebuffer for display. The simulator also avoids the CPU to GPU information transfer bottleneck by its storage of data on the video device memory. This is accomplished with the OpenGL Shading Language that is used through PyOpenGL to achieve a better visual description of information than is provided in other tools. A background on PyOpenGL and its comparison to other visualization libraries can be found in Rossant and Harris (2013).

Optimal Visualization Rate
The fourth approach to accelerate simulations is to display field matrices at a rate optimized for continuous visualization perception. Two images are perceived simultaneously when there is an interval of less than 30 ms between them (Wertheimer et al., 2012). The simulator takes advantage of the temporal lag in biological visual perception by stopping the numerical calculations to submit the field data to the GPU once within every 30 ms. This allows for the numerical part of the simulations to continue with fewer stoppages, resulting in faster simulations.

3D Visualization
The open source and cross-platform show3D library was written for the Neural Field Simulator to display field information. The library's visualization of neural fields expands two dimensional neural field data into a third dimension to better observe the differences in field locations. This is achieved by raising every

APPLICATIONS
The simulator can be used to analyze spatio-temporal neural field dynamics. The simulator's open source code allows modifications and extensions to be added to the code. The subsequent sections describe few of these possible applications.
Introducing finite axonal transmission speed in neural fields substantially slows numerical computation. However, to omit finite transmission speed is to neglect biological physiology (Idiart and Abbott, 1993). Hutt and Rougier (2010) have suggested to implement finite axonal transmission speed in a computationally efficient manner that is utilized by the simulator. Numerically, the speed is infinite if c≥l/ √ 2 t and there is increasing delay as c decreases.

Breather
Breather oscillations have been reported in theory (Folias and Bressloff, 2005;Hutt and Rougier, 2010) and experiments (Wang, 2010). As shown in Hutt and Rougier (2010), breathers are solutions of Equation (1) for finite axonal transmission speeds. They can be obtained and visualized in the simulator by assuming a temporally constant external input I in Equation (1). For a Gaussian-shape input with its apex at the center of the field, one overwrites the I variable section in the values.py file as sigma = 5.65685425 I = 20 * np.exp(-x ** 2/sigma ** 2)/(sigma ** 2 * np.pi) and change the showData variable assignment near the beginning of the values.py file to showData = 3 to show the input I in the simulation. In the definition of I, the space variable x ∈ R 2 is defined to cover the spatial domain (not shown in the code snippet). A field input similar to Figure 3A can be seen when the simulator is run.
An inhibitory synaptic connectivity kernel, K in Equation (1), can be implemented for the breather and viewed by changing the showData and K variables in the values.py file to showData = 4 K = -4 * np.exp(-x/3)/(18 * np.pi) and running the simulator. Here, x ∈ R 2 . An inhibitory synaptic kernel similar to the one in Figure 3B can be subsequently viewed.
After overwriting I and K as noted above, replace the following variables and function in the values.py file with the values below: after running the program.

Turing Patterns
Turing patterns (Turing, 1952) have been reported in neural field models Elvin et al., 2009; Steyn-Ross et al., With c = 6364 the effective speed is infinite for the given l and dt properties. Figure 5 shows the simulation starting with random field potential noise. A pattern begins to emerge at 1 s and evolves into a temporally constant Turing pattern after approximately 5 s.   were chosen for clear displays of different (vertical, cone, and horizontal) patterns.

Extensions
The simulator, being open source, allows the tailoring of code to provide modifications and extensions. An example extension is the addition of a graphical interface to modify parameters and simulate neural fields. Figure 8 shows an example interface coded with the wxPython 6 toolkit. Simulations are started, paused, continued, and started anew by clicking a button. Neural field parameters can be modified prior to and during simulations by clicking on the appropriate area of the interface and completing a pop-up dialogue. Running simulations are automatically paused when a mouse hovers above parameter areas of the interface. There is a symbiosis among the show3D library discussed in Section 2.3 and the parameter selection interface. It is possible to view the external input stimulus, kernel, and firing rate in the GLFW window by adding a mouse event and hovering over these sections to automatically view the respective matrices. Viewing these elements while changing their parameters can help to fine-tune field parameters. Moving the mouse away from these areas unpauses paused simulations 6 http://www.wxpython.org. Frontiers in Neuroinformatics | www.frontiersin.org and the field potential matrix is shown in the GLFW window. Further synergy between the interface and show3D library is implemented with the option to alter graph values from the interface where z axis limits and colors can be selected.

DISCUSSION
The Neural Field Simulator and its dependencies are crossplatform. However, the simulator interacts with graphics hardware using system-specific drivers which can result in problems on some operating systems. The graphical user interface in Section 3.4 is an example of this where the cross platform wxPython toolkit uses OpenGL to draw to the screen. The show3D library likewise uses OpenGL to interact with GPU. The graphical user interface and show3D library function symbiotically on Mac systems via the Apple Graphics Library. Conversely, on other operating systems such as Linux and Windows, unfortunately the separate utilization of OpenGL causes the simulator to crash. To this end, the current version of the simulator is released without the addition of extensions in order to operate properly on every major operating system. Nevertheless, a graphical interface can be a good choice with an appropriate single operating system.
Apart from the software implementation, in future work some model assumptions can be released. The square geometry can be recast easily to a rectangular geometry, whereas more general geometries (e.g., the impressive implementation work in The Virtual Brain Sanz Leon et al., 2015) will take more numerical effort. Periodic boundary conditions guarantee the simple application of the FFT, effective implementations of other boundary conditions like Dirichlet conditions [V(z, t) =const, z ∈ ∂ ] will demand some implementation changes in the spatial integral computation. Such modifications may still retain the fundamental implementation of the FFT. In contrast, rejecting the homogeneity hypothesis of spatial interactions, i.e., K = K(x, y) = K(x − y), abolishes the convolution structure and slows down the numerical simulation, cf. Appendix I.
Future work will extend the NFM model to multiple equations to render the model even more biologically plausible. In addition, an extension of the implementation to a mixture of constant and space-dependent delays as considered by Veltz and Faugeras (2011) will be interesting.