Simulations of pattern dynamics for reaction-diffusion systems via SIMULINK

Natural systems exhibit an amazing diversity of patterned structures in both living and nonliving systems[1], such as animal coats (e.g., zebra stripes, leopard spots and fish spirals), chemicals in a gel[2], laser light in a cavity[3], charges on the surface of a semiconductor[4], ecological balance between two species[5] and neuronal activations in human cortex[6]. Generally, there exist two kinds of natural patterns[7]:

Pattern dynamics research focuses on universal pattern-forming mechanisms for the latter case. Away from thermodynamic equilibrium, some experimentally observed 2D spatial patterns have been reported: Rayleigh-Bénard convection patterns[8]; reaction-diffusion Turing patterns[9]; Faraday-wave patterns[10]; vibratory granular patterns[11]; slime-mold spiral patterns[12] and Ca²⁺ spiral waves in Xenopus laevis eggs[13]. Although these patterns are different in the spatiotemporal scales or detailed pattern-forming mechanisms, they have similar structures.

More than half a century ago, British mathematician Alan Turing strove to understand the spontaneous processes generating these structures. In his famous 1952 paper[14], he proposed a reaction-diffusion model explaining potential mechanism for animal coats: At a certain stage of embryonic development, the reaction and diffusion between molecules, known as morphogens, and other reactors, lead to the breaking of symmetry of the homogeneous state. The morphogens spontaneously evolve to a non-uniform state, leading to the unique textures seen on animal skin. The key pattern-forming idea in Turing’s paper is the system’s spontaneous symmetry-breaking mechanism, also known as a Turing bifurcation (spatial instability).

However, for nearly 30 years Turing’s paper drew little attention for two reasons: morphogens had not been found in biological systems; negative chemical concentrations are permitted in Turing’s model, but are not accepted by chemists. From the late 1960s, a Brussels team led by Ilya Prigogine (winner of 1977 Nobel prize in chemistry) was endeavouring to prove Turing’s theory. Prigogine developed a reaction-diffusion model, Brusselator[15], to show the existence of Turing patterns that obey the rules of chemical kinetics. The Brusselator is a mathematical representation of the interaction between two morphogens, a reactor and an inhibitor, competitively reacting in time and spreading in space, which could give rise to a symmetry-breaking transition bifurcating from a homogeneous to patterned state, either stationary in a spatial (Turing) structure or in a temporal (Hopf) oscillation[15]. The Brusselator model further revealed the “secret” of Turing patterns: a coupling between nonlinear reaction kinetics and distinct diffusion rates such that the inhibitor should diffuse more rapidly than the activator[16].

This Brusselator proposed activator–inhibitor interaction is now a well known universal principle explaining regular pattern formation in chemical[17–19], ecological[5, 20] and physical[21, 22] systems, as well as biological systems[23, 24].

As pattern-forming systems essentially consist of coupled differential equations, the simulated patterns are time- and space-dependent solutions of the differential equations. The MATLAB built-in fourth-order Runge-Kutta solver ode45 and custom Euler methods are common ways to integrate differential systems. The implementation of these algorithms are well explained and demonstrated in many studies[25–28]. For example, in Yang’s book[25], at the end of Part II Yang presents a piece of concise MATLAB code for efficiently simulating simple reaction-diffusion systems. With some modifications, Yang’s programs can be used to simulate pattern formation in a wide range of applications of nonlinear reaction-diffusion equations. Some of these examples are discussed in detail in Part III of Yang’s book. Alternative to MATLAB, there are other options for pattern simulations such as MEREDYS[29], implemented in the Java programming language, interpreting the NEUROML model description language[30]. Besides, there are other programming environments applicable to modelling pattern formations such as COMSOL[31] and MODELICA[32].

Although it is efficient to solve differential equations in MATLAB or other programming platforms, their code script-based pattern simulations require high-level programming skills to tune the model parameters in order to examine the pattern dynamics; this limits their uptake by non-programmers. A few attempts for graphic-based pattern simulations have been made in the last decades. For example, READY[33] is an ideal program for getting started with reaction-diffusion systems. It runs fast (taking advantage of multi-core architectures), is easy to use (graphic-based) and runs on multiple platforms. In addition, there are other Java applets allowing easy pattern simulations such as Gray-Scott Simulator[34], showing the Gray-Scott pattern[35]: phenomena in grids of points that are connected “amorphously”. This closely models the development of a biological system of living cells. Similarly, Lidbeck has created another Java application[36] with extensive presets and options, and even allows 3-D of pattern simulations. However, these graphical user interface (GUI)-based softwares designed mainly for demonstrations of specific (pre-loaded) models. READY supports custom models but is written in READY-specific codes. So technically, READY is still a coding-based platform with a GUI interface. In the community of pattern dynamicists, there may be a demand for a software platform circumventing the programming burden required to implement numerical simulations of biologically-relevant pattern-forming systems.

The aim of this paper is to introduce SIMULINK modelling for simulating pattern dynamics. SIMULINK, an add-on product to MATLAB, provides an interactive, graphical environment for simulating and analysing dynamic systems. It enables modelling via a graphical programming language based on block diagrams. SIMULINK has been used extensively in engineering field[37] such as control theory and digital signal processing for multidomain simulation[38] as well as model-based design[39]. Besides, SIMULINK users have extended its applications in other areas such as medical device prototyping[40], heat transfer modelling[41, 42] and chemical reactions[43].

The present paper introduces the application of SIMULINK to pattern simulation studies, and it is organised as follows. We start, for pedagogical reasons, with a brief demonstration of the SIMULINK in modelling differential equations. Then, we construct the well-known Brusselator model in SIMULINK. By extending the Brusselator modelling strategies, we construct a SIMULINK-based human cortical model[44] (developed by the authors’ team) that incorporates the pattern-forming essentials while retaining neurophysiological accuracy: the cortical model comprises interacting populations of excitatory and inhibitory neurons which communicate via chemical (neurotransmitter-controlled) and electrical (gap-junction) synapses. In the model, the interaction between chemical kinetics and electrical diffusion allow for the emergence of a comprehensive range of patterns, which may be of direct relevance to clinically observed brain dynamics such as epileptic seizure EEG spikes[6], deep-sleep slow-wave oscillations[45] and cognitive gamma-waves[46, 47].

Finally, we examine the Brusselator and cortical model pattern dynamics. These simulation results are compared with the theoretical predictions.

The diffusion constant D _U,V [with units (length)²/time] is an important parameter indicative of the diffusion mobility. For a multi-component system, the higher the diffusivity, the faster the species diffuse into each other. Here, f _U,V (U,V) is typically a nonlinear function of concentrations U and V.

Solving a pattern-forming system in the form of Eq. (1) requires interpreting the differential operators for time ∂/∂ t and space ∇². In the following example, we first model the van der Pol oscillator in SIMULINK to explain how we interpret the differential operator on time.

Van der Pol oscillator

The van der Pol oscillator was originally developed by the Dutch electrical engineer and physicist Balthasar van der Pol[48]. The van der Pol oscillator was the first mathematical model proposed for the heartbeat, and it has also been used to simulate brain waves[49, 50]:

$\frac{d^{2}x}{d{t}^2}-\mu \left(1-x^2\right)\frac{\mathit{\text{dx}}}{\mathit{\text{dt}}}+x=0$

(2)

We wish to solve this equation for the case μ = 1 with initial conditions x(0) = 2 and d x/d t = 0 at t = 0. The traditional way to solve a second-order differential equation is to convert to a pair of coupled first-order differential equations:

$\begin{array}{l}\dot{x}=y\\ \dot{y}=\mu \left(1-x^2\right)y-x\end{array}$

(3)

We would now integrate these equations with time using the MATLAB numerical integrator ode45. This helps to form the link with the integration in SIMULINK.

We code the first-order van der Pol equations into a MATLAB function^a as follows:

To solve Eq. (3), we specify the coefficient μ, the initial conditions and the time-span over which the integration is to proceed; then pass these values, along with the name of the van der Pol function, to the Runge-Kutta solver ode45:

The calculated results are plotted in Figure1.

Alternatively, we may use the SIMULINK construction of Eq. (2), as shown in Figure2.

At a first glance, the interface for SIMULINK is completely different from the code sheet. In SIMULINK, all calculating elements are displayed by blocks. We select blocks from the SIMULINK library, then connect them to build a model.

The basic principle to model a differential equation in SIMULINK is to find the input and output of an integrator. Since we have:

$\int \left[\int \frac{d^{2}x}{d{t}^2}\mathit{\text{dt}}\right]\mathit{\text{dt}}=\int \frac{\mathit{\text{dx}}}{\mathit{\text{dt}}}\mathit{\text{dt}}=x$

(4)

then it follows that for a second-order differential equation, we need at least two integrators. As seen in Figure2, we first place an integrator block (the left block labelled with$\frac{1}{s}$) to process the inner integration of Eq. (4):$\int \frac{d^{2}x}{d{t}^2}\mathit{\text{dt}}$. The output of this integrator reads d x/d t, which is sent into the second integrator (the right block labelled with$\frac{1}{s}$). We assume the integrated x is known, thus being used to construct the input of the left integrator block, which is equivalent to$\frac{d^{2}x}{d{t}^2}$ with the form:

$\frac{d^{2}x}{d{t}^2}=\mu \left(1-x^2\right)\frac{\mathit{\text{dx}}}{\mathit{\text{dt}}}-x$

(5)

The product block (labelled with ×) in Figure2 combines (1 - x ²) and dx/dt. The result is amplified by a gain (triangle block, valued μ), then passed through a function block where x is subtracted. Here, the RHS of Eq. (5) is constructed.

Modelling a differential equation in SIMULINK requires forming a closed loop, where the integrated variables are fed back into the system. Evolution proceeds until reaching the desired final time. The scope block shows the real-time output of the two integrators; the scope can be placed anywhere to monitor the response of a sub-system. The Out1 and Out2 terminals send outputs of two integrators to the MATLAB workspace for further analysis. The results of this SIMULINK model are exactly the same as shown in Figure1.

Both MATLAB and SIMULINK allow fixed or self-adaptive (i.e., auto) time-steps for the Runge-Kutta solver^b. Figure3 shows that the discrepancy between MATLAB and SIMULINK Runge-Kutta solvers in either fixed or auto time-step mode are sufficiently small (< 10^-10). Consequently, we can see that the accuracy of the model simulation does not depend on the modelling platform since MATLAB and SIMULINK share the same integration algorithm to solve differential equations. However, modelling in SIMULINK is more straightforward and intuitive, and requires less programming skill than the MATLAB code sheet. The original mathematical equations can be converted into SIMULINK by matching its pattern with SIMULINK blocks directly. Moreover, in SIMULINK, by simply adding more blocks, or replacing blocks, a new model is able to be built in a very short time. SIMULINK may be an ideal tool to efficiently perform the simulations of a mathematical model. In the next section, we will extend the SIMULINK modelling method to describe a Brusselator system considering both its temporal and spatial evolutions on 1-D and 2-D Cartesian grids.

Readers should be aware of the choice of an appropriate differential solver for a specific problem, depending on the stiffness of differential equations. Applying a wrong solver may lead to either unstable solution or exceptional computation time. However, it is practically difficult to identify the stiffness of a differential model, thus one should try at least two different solvers, and compare the results. If they concur, i.e. give the same solution, they are likely to be correct. As suggested by MATLAB help file, it is worthwhile to try ode45 first since it is the most widely used method. For pattern-forming systems, we can also compare the numerical solution with the theoretical prediction to identify the applicability of the solver. For the demonstrated Brusselator and cortical models, ode45 and ode23 both work well and give rise to similar result; moreover, the numerical solutions match well with the theoretical predictions in emergent patterns (see Results and discussion section). So we choose ode45 solver to integrate the differential-equation models in this paper.

Brusselator model

The Brusselator model describes the competition of two chemical species in a chemical reaction, and is the simplest reaction-diffusion system capable of generating complex spatial patterns. The competition between two reactors and the introduction of diffusion satisfy the key requirements for pattern formation[14]. The pattern dynamics of the Brusselator has been comprehensively examined in de Wit[15] and other researchers’ work[51–54]. Here, our purpose is to introduce the SIMULINK modelling strategies.

SIMULINK version of Brusselator model

The simplest form of the model reads[55],

$\begin{array}{l}\frac{\partial }{\partial t}X=A-(B+1)X+X^{2}Y+D_X{\nabla }^{2}X\\ \frac{\partial }{\partial t}Y=\mathit{\text{BX}}-X^{2}Y+D_Y{\nabla }^{2}Y\end{array}$

(6)

where X and Y denote concentrations of activator and inhibitor respectively; D _X and D _Y are diffusion constants; A is a constant and B is a parameter that can be varied to result in a range of different patterns.

The LHS of Eq. (6) is a partial derivative on time since X and Y are functions of both time and space. At the RHS, the spatial derivative is represented by a Laplacian operator ∇². In the numerical simulation, the spatial dimension of the model is discretised into a grid by using the finite difference method. In the two-dimensional system the Laplacian with respect to the concentration field U in the node (i,j) is calculated along the x and y directions simultaneously:

${\nabla }^2{U}_{i,j}\approx \frac{{\mathrm{\Delta }}_x^2{U}_{i,j}}{h_x^2}+\frac{{\mathrm{\Delta }}_y^2{U}_{i,j}}{h_y^2}$

(7)

where

$\begin{array}{l}{\mathrm{\Delta }}_x^2{U}_{i,j}=U_{i+1,j}-2U_{i,j}+U_{i-1,j};\\ {\mathrm{\Delta }}_y^2{U}_{i,j}=U_{i,j+1}-2U_{i,j}+U_{i,j-1}\end{array}$

(8)

The h _x,y demonstrators in Eq. (7) are the respective x and y grid spacings; they define the spatial resolution. Assuming h ≡ h _x  = h _y (i.e., a square grid), the discrete Laplacian operation in a one-dimensional Cartesian coordinates along the y-axis has the form:

${\nabla }_{\text{1D}}^2{U}_{i,j}\approx \frac{U_{i,j+1}-2U_{i,j}+U_{i,j-1}}{h^2};$

(9)

for the two-dimensional case, we have

$\begin{array}{l}{\nabla }_{\text{2D}}^2{U}_{i,j}\approx \frac{U_{i+1,j}+U_{i-1,j}-4U_{i,j}+U_{i,j+1}+U_{i,j-1}}{h^2}\end{array}$

(10)

In SIMULINK, we initialise the Brusselator model as a column vector consisting of a 60 × 1 grid (spatial resolution = 1 cm/grid-point) for the one-dimensional case; or as a 60 × 60 grid for the two-dimensional case. Grid edges are joined to give toroidal boundaries.

The Laplacian operator${\nabla }_{\text{1D}}^2$ in Eq. (9) is implemented as a circular convolution of the 3 × 1 second-difference kernel$L_{\text{1D}}^y$ acting along the y-axis:

$\begin{array}{l}{\nabla }_{\text{1D}}^2\approx L_{\text{1D}}^y=\frac{1}{h_y^2}\left[\begin{array}{c}1\\ -2\\ 1\end{array}\right]\end{array}$

(11)

The two-dimensional Laplacian operator${\nabla }_{\text{2D}}^2$ in Eq. (10) is built up from the sum of two orthogonal L _1D operators:

$\begin{array}{l}{\nabla }_{\text{2D}}^2\approx L_{\text{2D}}=\frac{1}{h_x^2}\left[\begin{array}{ccc}0& 1& 0\\ 0& -2& 0\\ 0& 1& 0\end{array}\right]+\frac{1}{h_y^2}\left[\begin{array}{ccc}0& 0& 0\\ 1& -2& 1\\ 0& 0& 0\end{array}\right]=\frac{1}{h^2}\left[\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right]\end{array}$

(12)

where we have again assumed a square grid so that h _x  = h _y  = h.

In SIMULINK, the 1-D or 2-D Laplacian operator with toroidal boundaries is processed through two blocks: The “wrap-around” and “2-D CONV” (can process both 1-D and 2-D convolutions). The “wrap-around” block wraps the input matrix on both axes to allow a valid convolution in the “2-D CONV” block against the Laplacian kernel L to return the cyclic convolution. We created a subsystem to compute the convolution, as shown in Figure4.

In Figure4, the custom block labelled “wraparound” contains codes extracted from the convolve2() function^c.

The reason we introduce the custom block is that the SIMULINK built-in 2-D CONV block provides only the “valid” (non-flux) boundary condition, and cannot handle periodic boundaries.

Following the ideas of modelling the van der Pol oscillator, we can easily convert Eq. (6) to SIMULINK blocks, seen in Figure5.

SIMULINK versions of Waikato cortical model equations

The authors’ team at the University of Waikato (New Zealand) has developed a mean-field cortical model which encapsulates the essential neurophysiology of the cortex, while remaining computationally cost efficient[44]. The model envisions the cortex as a continuum in which pools of neurons are linked via chemical and electrical (gap-junction) synapses.

Model background

The Waikato cortical model has a rich history of development: The model is based on early work by Liley et al.[56], with enhancements motivated by Rennie et al.[57] and Robinson et al.[58]; and has been recently extended to include electrical synapses (also referred as gap junctions)[44, 47] supplementing neural communications via chemical synapses.

The Waikato cortical modelling assumptions are:

1. 1.

   Cortical element is the “macrocolumn” containing ∼100,000 neurons.

    
2. 2.

   There are only two distinct kinds of cortical neurons: 85% excitatory, and 15% inhibitory neurons. “Excitatory” and “inhibitory” are classified according to their effects on other neurons.

    
3. 3.

   There are three isotropic neuronal interactions: cortico-cortical (long-range from other macrocolumns), intra-cortical (short-range within macrocolumn) and connections from subcortical structures such as the thalamus and brain-stem. A macrocolumn with diameter ∼1 mm and depth ∼3 mm is sketched in Figure 6. We further assume that inhibitory connections via their short axons are locally restricted within a macrocolumn. In contrast, the axons from excitatory neurons are often longer and extensively branched, having length ranging from centimetres to several meters (e.g., in a giraffe’s neck), allowing exclusively excitatory cortico-cortical connections. Both excitatory and inhibitory connections are permitted in local neuron connections.

    
4. 4.

   Neurons exchange information via both chemical and electrical (gap-junction) synapses. Gap-junctions are more abundant within inhibitory populations, while being rare within excitatory neuronal populations.

    

The authors’ team first introduced gap-junctions into the cortical model in the paper Gap-junction mediate large-scale Turing structures in a mean-field cortex driven by subcortical noise[44]. The cortical model is expected to exhibit similar dynamics to a chemical reaction-diffusion system: The interaction between the excitatory and inhibitory neurons is analogous to to the competition between the activator and inhibitor of the chemical reaction-diffusion model, e.g. the Brusselator. The gap-junction strength between inhibitory neurons also plays the same role as the diffusion terms in the Brusselator, which allows a spatial evolution of the patterns. Consequently, it is reasonable that the cortical model exhibits similar pattern dynamics for those seen in the chemical reaction-diffusion system.

Model equations

The neuron populations deliver spike flux ϕ _ab from sources Q _a at a distance following damped wave equations[58]:

$\left[{(\partial /\partial t+v{\mathrm{\Lambda }}_{\mathit{\text{ab}}})}^2-v^2{\nabla }^2\right]{\varphi }_{\mathit{\text{ab}}}=v^2{\mathrm{\Lambda }}_{\mathit{\text{ab}}}^2{Q}_a$

(13)

The subscript ab is read “a → b”, the connection from the presynaptic neuron type a to postsynaptic neuron type b. a,b can be either e (excitatory) or i (inhibitory). Q _a is the spike-rate flux, which is a mapping from membrane voltage V _a to population-averaged firing rates:

$Q_a=\frac{Q_a^{\text{max}}}{1+e^{-C(V_a-{\theta }_a)/{\sigma }_a}}$

(14)

The total input flux arriving at the post-synapse is given as,

$M_{\mathit{\text{ab}}}=\underset{\text{long-range}}{\underbrace{N_{\mathit{\text{eb}}}^{\alpha }{\varphi }_{\mathit{\text{eb}}}^{\alpha }}}+\underset{\text{local}}{\underbrace{N_{\mathit{\text{ab}}}^{\beta }{\varphi }_{\mathit{\text{ab}}}^{\beta }}}+\underset{\text{subcortical}}{\underbrace{{\varphi }_{\mathit{\text{eb}}}^{\text{sc}}}}$

(15)

where superscripts α and β distinguish long-range and local chemical synaptic inputs;$N_{\mathit{\text{eb}}}^{\alpha }$ and$N_{\mathit{\text{ab}}}^{\beta }$ are the number of such input connections. The cortex is driven by subcortical noise which enters the intra-cortex. Subcortical excitation${\varphi }_{\mathit{\text{eb}}}^{\text{sc}}$ is modelled as small-amplitude spatiotemporal Gaussian-distributed white noise superimposing on a background excitatory “tone”$〈{\varphi }_{\mathit{\text{eb}}}^{\text{sc}}〉$ whose constant level is set via a subcortical drive scale-factor s:

${\varphi }_{\mathit{\text{eb}}}^{\mathit{\text{sc}}}(t)=s〈{\varphi }_{\mathit{\text{eb}}}^{\text{sc}}〉+{\xi }_{\mathit{\text{eb}}}(t)\sqrt{s〈{\varphi }_{\mathit{\text{eb}}}^{\text{sc}}〉}$

(16)

where ξ _eb is the Gaussian-distributed white-noise sources being delta-correlated in time and space. Inclusion of white-noise stimulation is motivated by physiological evidence that the cortex requires a continuous background “wash” of input noise to function normally.

The total input flux Φ _ab is the temporal convolution of the dendrite impulse response H(t) with the input flux M _ab :

$\begin{array}{ll}{\mathrm{\Phi }}_{\mathit{\text{ab}}}& =(\text{dendrite response})\otimes (\text{input flux})\\ =H_b(t)\otimes M_{\mathit{\text{ab}}}(t)\\ ={\int }_0^t{H}_b(t-t^{\prime })M_{\mathit{\text{ab}}}(t^{\prime })d{t}^{\prime }\end{array}$

(17)

where the dendrite dynamics are modelled by the alpha-function impulse response

$H_b(t)={\gamma }_b^{2}t{e}^{-{\gamma }_{b}t}$

(18)

Reducing γ _i is equivalent to increasing the time-to-peak (1/ γ _i ) for the hyperpolarising GABA (gamma-aminobutyric acid) response, as indicated in Figure7.

By taking derivatives, Eq. (17) becomes

${\left(\frac{\partial }{\partial t}+{\gamma }_a\right)}^2{\mathrm{\Phi }}_{\mathit{\text{ab}}}={\gamma }_a^2{M}_{\mathit{\text{ab}}}(t)$

(19)

Therefore, the flux received by target neuronal populations are:

${\left(\frac{\partial }{\partial t}+{\gamma }_e\right)}^2{\mathrm{\Phi }}_{\mathit{\text{eb}}}=\left[N_{\mathit{\text{eb}}}^{\alpha }{\varphi }_{\mathit{\text{eb}}}^{\alpha }+N_{\mathit{\text{eb}}}^{\beta }{Q}_e+{\varphi }_{\mathit{\text{eb}}}^{\mathit{\text{sc}}}\right]{\gamma }_e^2$

(20a)

${\left(\frac{\partial }{\partial t}+{\gamma }_i\right)}^2{\mathrm{\Phi }}_{\mathit{\text{ib}}}=N_{\mathit{\text{ib}}}^{\beta }{Q}_i{\gamma }_i^2$

(20b)

Finally, those incoming spikes perturb the soma resting potentials:

$\begin{array}{l}{\tau }_b\frac{\partial V_b}{\partial t}=V_b^{\text{rest}}-V_b+\underset{chemical\; synapses}{\underbrace{{\rho }_e{\psi }_{\mathit{\text{eb}}}{\mathrm{\Phi }}_{\mathit{\text{eb}}}+{\rho }_i{\psi }_{\mathit{\text{ib}}}{\mathrm{\Phi }}_{\mathit{\text{ib}}}}}+\underset{electrical\; synapses}{\underbrace{D_{1,2}{\nabla }^2{V}_b}}\end{array}$

(21)

We write D _bb as the diffusive-coupling strength between electrically adjoined e ↔ e, i ↔ i neuron pairs. To simplify the notation, we write (D ₁,D ₂) ≡ (D _ee ,D _ii ).

Symbol definitions for the cortical model are summarised in Table1.

Symbols	Description	Value	Unit
Λ e,b	Inverse-length scale for e→b axonal connection	4	cm-1
v	Axonal conduction speed	140	cm s-1
$Q_{e,i}^{\text{max}}$	Maximum firing rate	30, 60	s-1
θ e,i	Sigmoid threshold voltage	-58.5, -58.5	mV
σ e,i	Standard deviation for threshold	3, 5	mV
C	Constant	$\pi /\sqrt{3}$
$N_{\mathit{\text{eb}}}^{\alpha }$	Long-range e → b axonal connectivity	2000	-
$N_{\mathit{\text{eb}},\mathit{\text{ib}}}^{\beta }$	Local e → b, i → b axonal connectivity	800, 600	-
γ e	Excitatory rate-constant	170	s-1
γ i	Inhibitory rate-constant	50	s-1
$〈{\varphi }_{\mathit{\text{eb}}}^{\text{sc}}〉$	e→b tonic flux entering from subcortex	300	s-1
τ e,i	Neuron time constant	0.04, 0.04	s
$V_{e,i}^{\text{rev}}$	Reversal potential at dendrite	0, -70	mV
$V_{e,i}^{\text{rest}}$	Neuron resting potential	-64, -64	mV
$\mathrm{\Delta }{V}_{e,i}^{\text{rest}}$	Offset to resting potential	1.5, 0	mV
ρ e	Excitatory synaptic gain	1.00×10-3	mV s
ρ i	Inhibitory synaptic gain	-1.05×10-3	mV s
D 2	i ↔ i gap-junction diffusive coupling strength	0–2.0	cm2
D 1	e ↔ e gap-junction diffusive coupling strength	D 2/100	cm2

SIMULINK versions of Waikato cortical model equations

Let us first list the mathematical equations for the Waikato cortical model and examine their characteristics.

•  The cortico-cortical equation

  $\left[{\left(\frac{\partial }{\partial t}+v{\mathrm{\Lambda }}_{\mathit{\text{eb}}}\right)}^2-v^2{\nabla }^2\right]{\varphi }_{\mathit{\text{eb}}}^{\alpha }={(v{\mathrm{\Lambda }}_{\mathit{\text{eb}}})}^2{Q}_e$
  can be arranged by collecting temporal derivatives to the LHS:

  $\frac{{\partial }^2}{\partial t^2}{\varphi }_{\mathit{\text{eb}}}^{\alpha }+2v{\mathrm{\Lambda }}_{\mathit{\text{eb}}}\frac{\partial }{\partial t}{\varphi }_{\mathit{\text{eb}}}^{\alpha }=v^2{\nabla }^2{\varphi }_{\mathit{\text{eb}}}^{\alpha }-v^2{\mathrm{\Lambda }}^2{\varphi }_{\mathit{\text{eb}}}^{\alpha }+{(v{\mathrm{\Lambda }}_{\mathit{\text{eb}}})}^2{Q}_e$

  (22)
• The intra-cortical equations

  $\begin{array}{l}{\left(\frac{\partial }{\partial t}+{\gamma }_e\right)}^2{\mathrm{\Phi }}_{\mathit{\text{eb}}}=\left[N_{\mathit{\text{eb}}}^{\alpha }{\varphi }_{\mathit{\text{eb}}}^{\alpha }+N_{\mathit{\text{eb}}}^{\beta }{Q}_e+{\varphi }_{\mathit{\text{eb}}}^{\mathit{\text{sc}}}\right]{\gamma }_e^2\\ {\left(\frac{\partial }{\partial t}+{\gamma }_i\right)}^2{\mathrm{\Phi }}_{\mathit{\text{ib}}}=\left[N_{\mathit{\text{ib}}}^{\beta }{Q}_i\right]{\gamma }_i^2\end{array}$
  have different RHS, but their LHS are in the same mathematical pattern:

  ${\left(\frac{\partial }{\partial t}+\gamma \right)}^2\mathrm{\Phi }=\frac{{\partial }^2}{\partial t^2}\mathrm{\Phi }+2\gamma \frac{\partial }{\partial t}\mathrm{\Phi }+{\gamma }^2\mathrm{\Phi }$

  (23)
  We can move the term γ ²Φ to the RHS of the intra-cortical equations, then the LHS of the intra-cortical equations have the expression:

  $\frac{{\partial }^2}{\partial t^2}\mathrm{\Phi }+2\gamma \frac{\partial }{\partial t}\mathrm{\Phi }$

  (24)

  which is similar to the LHS of the cortico-cortical equation.

• The soma equation

  ${\tau }_b\frac{\partial V_b}{\partial t}=V_b^{\text{rest}}-V_b+({\rho }_e{\psi }_{\mathit{\text{eb}}}{\mathrm{\Phi }}_{\mathit{\text{eb}}}+{\rho }_i{\psi }_{\mathit{\text{ib}}}{\mathrm{\Phi }}_{\mathit{\text{ib}}})+D_{\mathit{\text{bb}}}{\nabla }^2{V}_b$
  can be re-arranged as

  $\begin{array}{l}\frac{\partial V_b}{\partial t}=\frac{1}{{\tau }_b}\left[V_b^{\text{rest}}-V_b+({\rho }_e{\psi }_{\mathit{\text{eb}}}{\mathrm{\Phi }}_{\mathit{\text{eb}}}+{\rho }_i{\psi }_{\mathit{\text{ib}}}{\mathrm{\Phi }}_{\mathit{\text{ib}}})+D_{\mathit{\text{bb}}}{\nabla }^2{V}_b\right]\end{array}$

  (25)

Following the ideas of SIMULINK modelling in van der Pol oscillator, we need two integrator blocks for Eqs. (22) and (24), and two convolution processing for Eqs. (22) and (25).

The strategy for modelling a large system is to focus on its subsystems first, then connect them together. The Waikato cortical model has three major parts: cortico-cortical, intra-cortical and soma equations. Figure8 shows how neuronal fluxes are transferred from one to another: cortico-cortical flux${\varphi }_{\mathit{\text{eb}}}^{\alpha }$ is delivered to the long-range targets Φ _ee and Φ _ei ; intra-cortical flux Φ _ee and Φ _ei , Φ _ie and Φ _ii merge into the soma equations. The output of the soma V _e is connected to source neurons to form the closed loop through the excitatory sigmoid function:

$Q_e=\frac{Q_e^{\text{max}}}{1+e^{-C(V_e-{\theta }_e)/{\sigma }_e}}$

In following sections, we detail the SIMULINK implementation of the three subsystems (drawn as three blocks in Figure8) of the cortical model.

Cortico-cortical flux

The SIMULINK based cortico-cortical block (see Figure9) is converted from Eq. (22). The flux-source Q _e is a mapping from the excitatory soma voltage sent via the Goto block, to the firing-rate received via the From block. After two integrations, signals will be passed to the excitatory synapses via port-1 (upper right corner) and inhibitory synapses via port-2 in the intra-cortex.

Intra-cortical flux

In SIMULINK modelling, we divide Eq. (16) into two parts: the constant (i.e. dc-level) excitatory background$s〈{\varphi }_{\mathit{\text{eb}}}^{\text{sc}}〉$ and the one-off kick$\sqrt{s〈{\varphi }_{\mathit{\text{eb}}}^{\text{sc}}〉}{\xi }_{\mathit{\text{eb}}}$. As illustrated in Figure10, we use a Clock block to count the iteration step. Once the counter is above one, the “switch” will turn off the kick, allowing only the constant excitation to enter the intra-cortex (removing the Clock block would allow on-going noise stimulus from the subcortex). The 2-D spatial white noise are generated by the Band-Limited White Noise block.

The intra-cortical model describes how post-synaptic fluxes evolve over time. In Figure11, the local fluxes (input via the From block) along with the long-range fluxes (from input-1 at the left, labelled as${\varphi }_e^{\alpha }$) and subcortical drive are summed, then filtered at the post-synaptic dendrite, thus forming the post-synaptic fluxes Φ _ee at the output port-1 (upper right corner). The Φ _ee and Φ _ei flux models have symmetric structure.

We assume that the cortico-cortical fibres are exclusively excitatory, thus there are no long-range inhibitory fluxes entering into the soma. Figure12 shows that the local inhibitory fluxes Φ _ie come from local source Q _i only. The Φ _ie and Φ _ii models also have symmetric structure.

Soma voltage

Figure13 presents the soma model of the excitatory neuronal group. The short-range fluxes are accumulated at the soma, from ports 1 and 2. The soma voltage V _e is converted to firing-rate Q _e locally in this sub-model (block labelled with Q _e sigmoid), then fed back into the cortico-cortical and intra-cortical models.

Finally, we connect all subsystems to form the completed Waikato cortical model, as illustrated in Figure14. It follows the flux flow-chart of Figure8, with the detailed SIMULINK block connections shown in Figure15. We argue that such model-based-design is an advantage for representing differential equations in SIMULINK. Although SIMULINK is useful for rapid prototyping, the SIMULINK implementation of the cortical model runs slower than our pre-coded Euler integration^d since it is time-consuming to interpret the MATLAB function wrapround (see Figure4) in SIMULINK. A 60×60 grid (side length 20 cm) 5-s cortical simulation takes ∼10 s via MATLAB Euler integration (fixed time-step 0.8 ×10^-3 s), while ∼40 s via SIMULINK (auto time-step and in accelerator mode). Thus, it is recommended to avoid using MATLAB functions in SIMULINK unless necessary.

The strategy for modelling differential equation models in SIMULINK is addressed in this paper. To construct a system of differential equations in SIMULINK, we can directly convert the mathematical terms and operators to the SIMULINK graphical block diagrams. The key idea for programming differential systems in SIMULINK is to form a closed loop, such that the solution of the system can evolve in this loop. The accuracy and reliability of the SIMULINK modelling method has been examined via comparing a van der Pol oscillator represented in MATLAB code-script and SIMULINK block-diagram, showing the SIMULINK model has comparable performance with the code-script version.

Using a well-known Brusselator model, we demonstrated two main SIMULINK modelling strategies for a reaction-diffusion system: interpretation of 2-D convolution with the periodic boundary condition; hybrid programming in SIMULINK with MATLAB functions. The pattern simulations of the Brusselator in SIMULINK are in good agreement with the predictions via bifurcation theories.

Following the pattern-forming theories, we introduced a cortical model with competitive neuronal interactions and diffusions. Unlike the simple Brusselator, the cortical model is a complicated system comprised of distinct cortical connections. Here, we showed how to build these subsystems in SIMULINK. Finally, we connected all subsystems to form the completed Waikato cortical model. The simulations of the cortical model exhibit Turing and mixed Turing–Hopf patterns, the clinical relevance of which is briefly discussed. As the main aim of this paper is to introduce SIMULINK modelling, readers are referred to Steyn-Ross et al.’s recent publications[6, 44–47] for comprehensive investigations of the cortical model.

^a vanderpoldemo is a MATLAB pre-coded function

^bfixed time-step ODE solvers are not built into MATLAB, but they can be acquired from a release by the MathWorks Support Team:http://www.mathworks.com/matlabcentral/answers/uploaded_files/5693/ODE_Solvers.zip

^cThe two-dimensional circular convolution algorithm was written by David Young, Department of Informatics, University of Sussex, UK. His convolve2() code can be downloaded from MathWorks File Exchange:http://www.mathworks.com/matlabcentral/fileexchange/22619-fast-2-d-convolution

^d MATLAB simulation codes were written by Alistair Steyn-Ross. The complete codes, plus README files and movies of cortical dynamics, are available from the web site:http://www2.phys.waikato.ac.nz/asr/