Open-source tools for dynamical analysis of Liley's mean-field cortex model

Mean-field models of the mammalian cortex treat this part of the brain as a two-dimensional excitable medium. The electrical potentials, generated by the excitatory and inhibitory neuron populations, are described by nonlinear, coupled, partial differential equations, that are known to generate complicated spatio-temporal behaviour. We focus on the model by Liley {\sl et al.} (Network: Comput. Neural Syst. (2002) 13, 67-113). Several reductions of this model have been studied in detail, but a direct analysis of its spatio-temporal dynamics has, to the best of our knowledge, never been attempted before. Here, we describe the implementation of implicit time-stepping of the model and the tangent linear model, and solving for equilibria and time-periodic solutions, using the open-source library PETSc. By using domain decomposition for parallelization, and iterative solving of linear problems, the code is capable of parsing some dynamics of a macroscopic slice of cortical tissue with a sub-millimetre resolution.


Introduction
Models of cortical dynamics come in two main families: network models and mean-field models.The former describe many interacting neurons, each with their own dynamical rules, while the latter describe electrical potentials, generated collectively by many neurons, as continuous in space and time.These potentials can be thought of as averages over a number of macrocolumns, groups of hundreds of thousands of neurons in columnar structures at the surface of the cortex.The reason to abandon the description of individual neurons and pass to the mean-field limit, in analogy to the thermodynamic limit in statistical physics, is twofold.Firstly, the description and analysis of a substantial piece of the cortex with a network model is not tractable since it would contain billions of neurons, and many times more connections between them.As demonstrated by recent publications, such as by Izhikevich and Edelman [1] or by the Blue Brain team [2], progress in super computing allows for the simulation of ever larger neuronal networks, that reflect actual brain dynamics.However, it is hard to see how the output of such models can be analysed, other than by purely statistical techniques.In contrast, mean-field models can be analysed as infinite-dimensional dynamical systems.
The second advantage of the mean-field approach is that the electrical potentials, which appears as dependent variables, are directly related to the electroencephalograph (EEG) [3].The EEG is usually measured with electrodes on the scalp or, in exceptional circumstances, directly on the surface of the brain.In either case, the measured signal is not that of individual neurons, but that of many neurons, spread out over a few square centimetres or millimetres.Thus, the way the signals of individual neurons are smeared out by the spatial averaging of mean-field modelling is similar to the way they are mixed up in EEG measurements.Because of the direct link between the local mean potential and the EEG, mean-field models are sometimes called EEG models.
The origin of mean-field modelling lies in the nineteen seventies, when pioneers like Walter Freeman [4], Wilson & Cowan [5] and Lopes da Silva [6] started to model components of the human cortex with continuous fields.Over the past four decades, mean-field models have been used to study a range of open questions in neuroscience, such as the generation of the alpha rhythm, 8-13 Hz oscillations in the EEG (see, e.g., [6,7]), epilepsy (see, e.g., [8,9,10]) and anaesthesia [11].Also, they are used in models for sensorimotor control, pattern discrimination and target tracking [12].
Although mean-field models have been used in all these contexts, little analysis has been done on their behaviour as spatially extended dynamical systems.In part, this is due to their staggering complexity.The Liley model [13] considered here, for instance, consists of fourteen coupled Partial Differential Equations (PDEs) with strong nonlinearities, imposed by coupling between the mean membrane potentials and the mean synaptic inputs.The model can be reduced to a system of Ordinary Differential Equations (ODEs) by considering only spatially homogeneous solutions, and the resulting system has been examined in detail using numerical bifurcation analysis (see [14] and references therein).In order to compute equilibria, periodic orbits and such objects for the PDE model, we need a flexible, stable simulation code for the model and it linearization that can run in parallel to scale up to a domain size of about 2500 cm 2 , the size of a full-grown human cortex.We also need efficient, iterative solvers for linear problems with large, sparse matrices.In this paper, we will show that all this can be accomplished in the open-source software package PETSc [15].Our implementation consists of a number of functions in C that will be available publicly [16].
The goal of this computational work is similar to that of Coombes et al., who analysed "spots": rotationally symmetric, localized solutions in a model of a single neuron population in two dimensions [17].The study of such special solutions will help us parse the spatio-temporal dynamics of mean-field models.We will attempt to compute periodic orbits and other special solutions in a full-fledged, two-population mean field model without imposing any spatial symmetries.

Liley's model
The model we use was first proposed by Liley et al. in 2002 [13].The dependent variables are the mean inhibitory and excitatory membrane potential, h i and h e , the four mean synaptic inputs, originating from either population and connecting to either, I ee , I ei , I ie and I ii , and the excitatory axonal activity in long-range fibers, connecting to either population, φ ee and φ ei .The model equations are where index k = {e, i} denotes excitatory or inhibitory.The meaning of the parameters, along with some physiological bounds and the values used in our tests, are given in Table 1.A detailed description of these equations can be found in references [13,14].Here, we will focus on the aspects of the model most relevant for the numerical implementation.
There are two sources of nonlinearity, related to the coupling of the synaptic inputs to the membrane potential and vice versa.The former connection is quadratically nonlinear, while the latter is given by the sigmoidal function S, which describes the onset of firing as the potential exceeds the threshold value µ i,e .These nonlinearities tend to form sharp transitions of the potentials across the domain.That is one reason why we opted for a finitedifference discretization over a pseudo-spectral approach.Spectral accuracy would be of limited value in the presence of steep gradients and the finitedifference scheme can be parallelized much more efficiently.The second reason is that we would like to be able to change the geometry of the domain and the boundary conditions in future work.The finite-difference scheme is more flexible in that respect.
The only spatial derivatives in the model are those in the equations for the long-range connections.These are damped wave equations.We will discretize the Laplacian using a five-point stencil on a rectangular grid.In previous work, Bojak & Liley chose a second-order centered difference scheme for the time derivatives [11].A disadvantage of this approach is that the stability condition of this scheme dictates that we set the time step inversely proportional to the grid spacing.In practice, they used a time step of 0.05 ms.To avoid this obstacle, we implemented the unconditionally stable implicit Euler method, as described in Sec. 3.
Other authors have used this model with an additional diffusive term in the equations for the membrane potentials to model gap junctions [18].Inclusion of these terms can drastically change the bifurcation behaviour, as they can cause Turing transitions to space-dependent equilibria.Without the additional terms, a Hopf bifurcation from a spatially homogeneous equilibrium to a space dependent periodic orbit or a saddle-node bifurcation of this equilibrium can be the primary instability.The gap junction terms can readily be included in our implementation, and in Sec. 5 we will describe how to solve for equilibrium states that may depend on space.
We will test our implementation by comparing to, and extending, the computations of oscillations with a 40 Hz component by Bojak & Liley [19].The corresponding parameter values are listed in Table 1.The 40 Hz os-cillations arise spontaneously if the number of local inhibitory-to-inhibitory connections is changed slightly.We introduce a scaling parameter r by replacing N β ii → rN β ii .This is the only parameter that will be varied in our tests.

PETSc overview
Rather than creating our code from scratch, we opted to work with the Portable, Extensible Toolkit for Scientific Computation (PETSc): an opensource, object oriented library that is designed for the scalable solution and analysis of PDEs [20,15].PETSc is written in the C language, and is usable from C/C++ as well as Fortran and Python.We use PETSc in conjunction with the Scalable Library for Eigenvalue Problem Computations (SLEPc) [21], for the computation of eigenspectra of equilibrium and periodic solutions.Since our implementation uses some features of PETSc and SLEPc that are recent additions and are still being tested, we use the development version of both projects.
PETSc is split up into multiple components to address the various problems associated with solving PDEs numerically.For our purposes, we treat the DM component, which handles the topology of the discretization, as the most fundamental, from which we can easily derive memory allocation and communication for distributed vectors (Vec) and matrices (Mat).With vectors and matrices, we can now solve linear systems, such as those that arise in Newton iteration for implicit time-stepping and the computation of equilibria and periodic orbits.PETSc's component for this is called KSP, and it has numerous iterative solvers implemented, as well as preconditioners, (PC), to increase convergence rates.For implicit time-stepping, for example, we use GMRES , preconditioned with incomplete LU (ILU) factorization, combined with the block Jacobi method [22,23].On top of the linear solvers come the nonlinear solvers, PETSc's SNES component, which implements a few different methods, such as globally convergent Newton iteration with line search [24].Finally, PETSc provides a timestepping component, TS, to obtain time dependant solutions.Implemented here are numerous explicit and implicit schemes such as adaptive stepsize Runge-Kutta and implicit Euler.The implicit schemes use the KSP component.A schematic of the hierarchy discussed here can be found in Fig. 1.
For our dynamical systems calculations we will frequently need to compute specific eigenvalues and eigenvectors for system-sized matrices.For this end, we use SLEPc, which implements iterative eigenvalue solvers using PETSc Vec and Mat distributed data structures.The component of SLEPc that we use is EPS, which has a few algorithms for iteratively solving eigenproblems.Its default algorithm is Krylov-Schur iteration.

Geometry
Following earlier work by Bojak & Liley (e.g.[11,19]) we consider the PDEs on a rectangular domain with periodic boundary conditions.On this domain, we use a rectangular grid of N x by N y points.In the tests presented in Sec. 7, the domain and the grid are square.PETSc allows for more complicated domain shapes and grids, that can be encoded in the DM component, independent of the higher-level components.
Within DM, PETSc provides a simpler subcomponent, DMDA, for working with finite differences on structured grids such as our rectangle.If we specify a stencil to use for the spatial derivatives in the DMDA, PETSc will automatically handle numerous things for parallel execution, such as memory allocation and the communication setup for distributed vectors and for the distributed Jacobian matrix.

Fields
To make use of PETSc's solvers, the model must be written as a system of equations that is first order in time.This we achieve by introducing new states J jk and ψ ek according to with indices j, k = {e, i}.
We opted to use a struct, seen in Code.2.1, to store the fields, rather than a triply indexed array.This allows the code to be more readable in the function and Jacobian evaluation routines.For example, one accesses the φ ee component at grid point (x i , y j ) simply as u[j][i].phiee, provided that the elements of the array (Field **u;) are stored on the processor in which the call is made.

Parameters
All of the model parameters are stored in a struct designated as the application context.The application context is how PETSc gets problem related parameters into the user-defined functions needed by its solvers.Similar to the fields, this allows readable code for the parameters.For example, Code 2.1 Struct for the fields.one accesses the Γ ie parameter as user->Gamma ie, if user is defined as the pointer AppCtx *user;.How the parameters show up in our struct for the application context is shown in Code 2.2.

User supplied functions
In addition to the structs listed above, we need to provide PETSc with (at least) a C function that computes the vector field for a given state.We call this function FormFunction, and from this PETSc is capable of approximating the Jacobian with various finite difference methods.However, we also supply a C function to explicitly compute the Jacobian, named FormJacobian, because this allows for more efficient calculations, especially when looking at stepping the variational equations in Sec. 4.

Timestepping
We use the implicit Euler method to time-step the discretized equations.As mentioned in Sec.1.1, this allows us to take larger time steps than feasible with explicit methods.Since we are aiming to compute periodic orbits, rather than to generate long time series, the first order accuracy of the method is not an issue.Once a periodic orbit is computed, the time-step size can be reduced to increase accuracy.Another option is to use Richardson extrapolation to increase the order of accuracy, using the same nonlinear solving as described below.

Mathematical basis
We symbolically write the dynamical system as where N is the total number of unknowns after discretization, in our case 14 × N x × N y .The implicit Euler scheme for time integration is given by where the subscript represents the step number, dt the step size, and u 0 the initial conditions.This nonlinear equation is solved by Newton iteration: where the superscript denotes the Newton iterate, and du k is the solution to the linear system where ∂f /∂u denotes the N × N Jacobian matrix.Provided that the initial approximation, u 0 n+1 , is close enough to the actual solution of equation (11), this iteration should converge quadratically.This is achieved by making the initial approximation the result of an explicit Euler step As we scale up the size of our problems, it becomes the linear solve in equation ( 13) that takes most time.This problem is handled by using GMRES to solve the linear system.For large time steps, the spectrum of the matrix in Eq. 13 is spread out, and we need to precondition it for iterative solving.We make use ILU, which has shown to be reliable [25,26] for this type of problem.If we use more than one processor, PETSc uses distributed storage for the matrix, and combines ILU with block Jacobi preconditioning.
Code 3.1 PETSc code for setting up and running the timestepping.FormFunction and FormJacobian are user provided functions that compute the rhs of equation (10), and its Jacobian respectively.J is an appropriately allocated matrix to hold the Jacobian, and u a vector to hold the solutions.

Implementation
PETSc provides a simple interface for timestepping in its TS component.The basic code required to set up a TS is given in Code 3.1.With a TS set up like this, the timestepping parameters are set from command line arguments at run time.For example, to do implicit Euler timestepping for 40.67 ms with a time step of 0.1 ms, one needs to provide the arguments -ts type beuler -ts dt 0.1 -ts final time 40.67.In this specific case, since the final time is not an integer number of timesteps, PETSc will step past it, and interpolate at the desired time.

Mathematical basis
The variational equations for the dynamical system are written as and must be integrated simultaneously with the dynamical system (10).Solving the variational equations allow us to compute the stability of solutions, and is also an essential ingredient for the treatment of boundary value problems such as those that arise in the computation of periodic orbits.
Performing implicit Euler timestepping on the variational equations ( 15) requires solutions of the linear problems Since we already have the Jacobian of the dynamical system at timestep n + 1, stepping the variational equations requires only one additional N × N linear solve per time step.

Implementation
In PETSc, we implement the timestepping of the variational equations as a MATSHELL.A MATSHELL allows users to define their own matrix type.Within a MATSHELL, one needs to give a context for storing the relevant data and write functions for the desired matrix operation.For example, we point the operation MATOP MULT to a function that takes the initial state of the variational system v(0) as input, and outputs the result v(T ) at the end of the timestepping.The context we use for the time stepping of the variational equations is shown in Code 4.1.The function we provide for MATOP MULT works by first taking a step of the TS, then loading the Jacobian computed from that step and solving equation ( 16).This is repeated until the TS reaches its end.The MATSHELL thus defined can be used by SLEPc for the iterative computation of eigen pairs.In particular, we will use this approach to compute the Floquet multipliers of periodic orbits.

Equilibria
Having set up the function FormFunction for the right hand side of the dynamical system, and its Jacobian computation FormJacobian, also used for time integration, we can set up equilibrium calculations using PETSc's SNES component with very little effort.

Mathematical Basis
Equilibrium solutions to the dynamical system (10) are solutions that satisfy To solve this, we can set up a Newton iteration scheme with du coming from the solution of the linear system As with the timestepping, if the initial guess is good enough this will converge quadratically provided that ∂f ∂u u k is nonsingular.Unlike the case of time stepping, though, we do not always have a way to produce an initial approximation that is good enough.For stable equilibrium solutions, we can use timestepping to get close to an equilibrium, but this will not work for unstable equilibria.One possible solution is using globally convergent Newton methods.Using such methods we can find equilibria from very coarse initial data, at the cost of computing many iterations.The line search algorithm and the trust region approach (see, e.g.[24]) are implemented in the SNES component.
Stability of equilibrium solutions follows from the spectrum of the Jacobian.Because of the spatial symmetries of the model, these will mostly appear in groups.On a square domain, for instance, a single eigenvalue will be associated with up to eight eigenvectors, with wavenumbers (±k x , ±k y ) and (±k y , ±k x ).

Implementation
Setting up and using a nonlinear solver within PETSc is straightforward, as shown in Code 5.1.The default algorithm used by SNES is Newton's method with line search.
Code 5.1 Code snippet for solving for equilibria.Vectors r and u are preallocated, with u being the initial approximation, and J a preallocated matrix for the Jacobian.

Periodic solutions
The primary instability in the Liley model is often a Hopf bifurcation, and periodic orbits have been shown to play an important role in the dynamics of ODE reductions of the model (e.g.[14,27]).However, space dependent periodic orbits have not previously been computed and studied.Using PETSc data structures for bordered matrices, in conjunction with a MATSHELL structure, we can solve for periodic orbits based on the time stepping described in Secs. 3 and 4.

Mathematical basis
Periodic orbits solve the boundary value problem where φ is the flow of the dynamical system (10), and T is the period.Our strategy for solving this equation is essentially that of Sanchez et al. [28], namely Newton iterations combined with unconditioned GMRES iteration.Linearising Eq. 20 gives where D u φ is a matrix of derivatives of the flow with respect to its initial condition.Upon convergence, this is the monodromy matrix of the periodic orbit.This results in N equations in N + 1 unknowns, which must be closed by a phase condition.We opted for the use of a Poincaré plane involving one of the state variables, u k : where C is set appropriately, for instance to the time-mean value of u k .This choice gives the following bordered system where [D u φ(u, T )] k,.denotes the k th row of the matrix D u φ.An update can then be made to the approximate solution as The matrix D u φ is dense, so we should avoid calculating and storing it explicitly.Iterative solving of the linear problem, (23), requires the computation of matrix-vector products, which are constructed from the integration of the variational equation ( 15) with v = du and the vector field f (φ(u, T )) at the end point of the approximately periodic orbit.Since the governing PDE is dissipative, most of the eigenvalues of the monodromy matrix are clustered around zero.This aids the convergence of GMRES, without any preconditioning.Sanchez et al. [28] provide bounds for the number of GMRES iterations for the Navier-Stokes equation, and the convergence we observe for the Liley model is qualitatively similar.

Implementation
The problem of creating a bordered matrix system in a distributed environment is not a trivial one.The specific case that we have is one vector, u, that is sparsely connected and distributed among processors, and one parameter, T , that must exist and be synchronized across all processors.
PETSc's DM module has some recently introduced functionality that allows us to handle this in a straightforward way, letting us make use of the DMDA already used in the other types of calculations.
DMRedundant can be used for the T component of our extended system, as it has the precise behaviour that we require.Next, we use a DMComposite to join together the DMDA of the grid with the DMRedundant of the period.We can then derive vectors from this DMComposite, and use these vectors for PETSc's iterative linear solvers.PETSc code that illustrates this idea is shown in Code 6.1.
The matrix multiplication is done through a MATSHELL, and the struct that holds the relevant data is found in Code 6.2.Code 6.1 Additional DM pieces for extended vectors as in equation ( 24), assuming that da is the DM associated with the grid structure.The numerical arguments in DMRedundantCreate represent the processor where the redundant entries live (in global vectors), and the number of redundant entries respectively.

Example calculations
In this section, we present some computations that serve to validate our implementation and to investigate its efficiency.All tests are based on the parameter set in Table 1, and the scaling of the number of local inhibitoryto-inhibitory connections, r, is varied around the first bifurcation from an equilibrium to more complicated, spatio-temporal behaviour.
Fig. 2 shows the neutral stability curve for the spatially homogeneous equilibrium, which is the unique attractor of the model at small values of r.The primary transition is a Hopf bifurcation with spatial wave numbers that depend on the system size.For systems smaller than 2 × 2 cm 2 , the emerging periodic orbit is spatially homogeneous.For larger systems, space dependent orbits emerge, and their typical length scale converges to about 9.3 cm for large system sizes.These stability curves were computed by solving small eigenvalue problems for each combination of wavenumbers, independent from the PETSc implementation.The eigenvalues computed by Krylov-Schur iteration in SLEPc, presented in Sec.7.2, are in good agreement.
A partial bifurcation diagram, for spatially homogeneous solutions only, is shown in Fig. 3.In this diagram, the Hopf bifurcation is subcritical, and time series analysis indicates that the Hopf bifurcations associated with nonzero wave numbers are, too.The time series presented in Sec.7.1 was generated by starting from the equilibrium at r = 1 and adding a finite-size perturbation in the least stable direction, with wave number Figure 2: Neutral stability curve for the spatially homogeneous equilibrium of the Liley model with parameters set according to Table 1.Shown is the scaling parameter, r, versus the linear domain size, L, and wave numbers k = (k x , k y ) are shown in parenthesis.When varying r, only for very small domains the primary instability is spatially homogeneous.For domain sizes over 12.5 × 12.5cm 2 the location of the primary instability approaches r = 1.04 and the length scale of the leading instability approaches L/ k = 9.3 cm.

Timestepping
For the timestepping demonstration, we used a system size of 12.8 × 12.8 cm 2 with 0.5 mm resolution, resulting in a 256 × 256 grid, and N = 917, 504 unknowns in total.Setting the parameter r = 1.0, we initialize with the stable equilibrium solution perturbed by its least stable eigenmode, shown in Fig. 7. Since the equilibrium solution is stable, small perturbations just decay, but sufficiently large perturbations grow.The snapshots of Fig. 4  show behaviour that is nearly periodic, with a dominant period of 40Hz, as demonstrated by the power spectrum shown in the last panel.
Since the time-stepping code lies at the core of the periodic orbit solver, we carefully investigated its scaling with an increasing number of processors.Doubling the domain size, while keeping the grid spacing fixed, gives a dynamical system with N = 3, 670, 016 degrees of freedom.We time-stepped this system on a small subcluster of 2.4 GHz AMD Opteron nodes with gigabit interconnects.Apart from some minor load-balancing effects, the scaling is linear up to 16 processors, despite the relatively slow interconnects.

Equilibrium
We computed the whole equilibrium curve of Fig. 3 through parameter continuation, which is a trivial extension of the algorithm for computing equilibria, presented in Sec. 5.For each computed equilibrium solution, we took the Jacobian and used SLEPc to compute the eigenvalues with the largest real parts.The result is shown in Fig. 7.As predicted by the neutral stability curve computation, the (1, 1) mode turns unstable first, immediately followed by the (1, 0) mode.Around r = 1.08, the (0, 2) mode crosses the (0, 1) mode and proceeds to become the most unstable mode for larger values of r.The least stable eigenmode for r = 1.046, just after its eigenvalue has crossed zero, is shown in Fig. 7.

Periodic solutions
We tested the computation of periodic orbits on a smaller grid, namely 16 × 16 points, still with 0.5 mm resolution, and with r = 1.2.The primary Hopf bifurcation is sub critical, so there is no easy way to compute the branch of space-dependent periodic solutions.Instead, we computed one  of the spatially homogeneous orbits, for which an approximate solution can readily be obtained from analysis of the ODE reduction of the model.In fact, the upper part of the branch of periodic orbits shown in Fig. 3 is stable to all spatially homogeneous perturbations.
Starting from a coarse initial approximation, the Newton iterations converged faster than linear, and each Newton step took between 8 and 11 GMRES iterations, out of a maximum of N + 1 = 3585.Subsequently, we computed the most unstable multipliers, using SLEPc with the MATSHELL that computes products with the bordered matrix, as described in Sec.6.2.The most unstable multiplier is µ 1 = 1.111 and corresponds to a wave number (1, 1) perturbation.

Conclusion and future improvements
In the current paper, we have presented the basic implementation of the model and example computations to validate it and test its performance.The code will be available publicly [16].As it is built on top of PETSc, the user has access to a range of nonlinear and linear solvers and preconditioners, which can be used to solve the boundary value problems that typically arise in dynamical systems analysis.The periodic orbit computation, presented in Sec.7.3, is a simple example of such a boundary value problem, that has all the ingredients: a module for time-stepping the system and perturbations and a representation of user-specified, bordered matrices.
The next step in the development of the code is the implementation of pseudo-arclength continuation of equilibria and periodic orbits.This will   enable us, for instance, to complement the bifurcation diagram of the current test case, Fig. 3, with the branches of space-dependent periodic solutions that actually regulate the observed dynamics, in contrast to the highly unstable spatially homogeneous periodic orbits computed from an ODE reduction of the model.We expect that our implementation will be useful to researchers studying the dynamics of the Liley model, or similar models, such as the model with gap junctions proposed by Steyn-Ross et al. [18].Also, it could be useful for those who incorporate a similar mean-field model in, for instance, the control of robotic motion or network models of brain activity.

Figure 1 :
Figure 1: Schematic representation of the components of PETSc and SLEPc used in our code, and their relative hierarchy.

Code 4 . 1
The MATSHELL context for timestepping of the variational equations.The TS holds the relevant info for stepping the dynamical system

20 %Figure 3 :
Figure3: Partial bifurcation diagram showing the primary transition from a spatially homogeneous equilibrium to a space and time dependent periodic orbit.On the vertical axis the scaling parameter r is plotted, and on the vertical axis the (maximum of) the excitatory membrane potential.The branch of periodic solutions shown with a dashed line is a spatially homogeneous branch that is unstable to space-dependent perturbations.

Figure 4 :
Figure 4: Three snapshots of the excitatory membrane potential, 6 ms apart, of a solution at r = 1, near the primary Hopf bifurcation.The domain size is 12.8 × 12.8 cm 2 , the resolution is 0.5 mm and the time-step size 1 ms.The fourth panel shows the power spectrum of h e , averaged over the region inside the black square.
e a r s c a l i n g

Figure 5 :
Figure 5: Wall time for the computation of 50 time steps of 1 ms each on a 25.6 × 25.6 cm 2 domain with 0.5 mm resolution.The fully implicit Euler steps are computed with Newton iterations, each of which is solved for by GMRES, preconditioned with a combination of block Jacobi and ILU.The initial guess is given by an explicit Euler step.Two or three Newton iterations are sufficient to reduce the residual by a factor of 10 8 .About 80 Krylov vectors are computed by GMRES to bring the relative residual down to 10 −5 .The number of unknowns is N = 3, 670, 016.

Figure 6 :
Figure 6: The real parts of the leading two eigenvalue pairs of the spatially homogeneous equilibrium tracked in the scaling parameter r, for system size L = 12.8 cm.The primary transition is tied to wave numbers |k x | = |k y | = 1.The other curves shown are for wave numbers k x = 0, k y = ±1 and k x = ±1, k y = 0 and for k x = 0, k y = ±2 and k x = ±2, k y = 0.

Figure 7 :
Figure 7: The real part of the least stable eigenmode of the stable equilibrium located at r = 1.046.Displayed are the excitatory (left) and inhibitory (right) membrane potentials.The eigenvector, with wave numbers (1, 1), was computed by Arnoldi iteration and is scaled to have unit L 2 norm.

Figure 8 :
Figure 8: Residuals of the Newton iteration (left) and the corresponding GMRES iterations (right).The latter is normalised by the norm of the right hand side of Eq. 23, i.e. the Newton residual.The tolerance was set at 10 −8 for the Newton iteration and to 10 −5 for the GMRES iteration.Note the super linear convergence of the former.
Code 2.2 Application context struct with the model parameters.

Table 1 :
[19]ing, ranges and values for the model parameters.The values used for the tests presented in Sec.7 are taken from reference[19].