Simulations of pattern dynamics for reactiondiffusion systems via SIMULINK
 Kaier Wang^{1},
 Moira L SteynRoss^{1},
 D Alistair SteynRoss^{1}Email author,
 Marcus T Wilson^{1},
 Jamie W Sleigh^{2} and
 Yoichi Shiraishi^{3}
https://doi.org/10.1186/17520509845
© Wang et al.; licensee BioMed Central Ltd. 2014
Received: 19 February 2014
Accepted: 7 April 2014
Published: 11 April 2014
Abstract
Background
Investigation of the nonlinear pattern dynamics of a reactiondiffusion system almost always requires numerical solution of the system’s set of defining differential equations. Traditionally, this would be done by selecting an appropriate differential equation solver from a library of such solvers, then writing computer codes (in a programming language such as C or Matlab) to access the selected solver and display the integrated results as a function of space and time. This “codebased” approach is flexible and powerful, but requires a certain level of programming sophistication. A modern alternative is to use a graphical programming interface such as Simulink to construct a dataflow diagram by assembling and linking appropriate code blocks drawn from a library. The result is a visual representation of the interrelationships between the state variables whose output can be made completely equivalent to the codebased solution.
Results
As a tutorial introduction, we first demonstrate application of the Simulink dataflow technique to the classical van der Pol nonlinear oscillator, and compare Matlab and Simulink coding approaches to solving the van der Pol ordinary differential equations. We then show how to introduce space (in one and two dimensions) by solving numerically the partial differential equations for two different reactiondiffusion systems: the wellknown Brusselator chemical reactor, and a continuum model for a twodimensional sheet of human cortex whose neurons are linked by both chemical and electrical (diffusive) synapses. We compare the relative performances of the Matlab and Simulink implementations.
Conclusions
The pattern simulations by Simulink are in good agreement with theoretical predictions. Compared with traditional coding approaches, the Simulink blockdiagram paradigm reduces the time and programming burden required to implement a solution for reactiondiffusion systems of equations. Construction of the blockdiagram does not require highlevel programming skills, and the graphical interface lends itself to easy modification and use by nonexperts.
Keywords
Background
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]:

Thermodynamic equilibriumbased, such as crystal structures in inorganic chemistry and spontaneously emergent organic polymer patterns. The forming mechanisms of these patterns have been extensively studied and wellexplained via thermodynamics and statistical physics: When a system’s temperature decreases, its entropy and the Gibbs free energy accordingly become smaller (systems tend to move towards a thermodynamic equilibrium where the Gibbs free energy reaches a minimum). Thus the principle of minimum energy allows the system to form certain spatial structures.

Thermodynamics farfromequilibrium, such as examples mentioned at the beginning of this section. These patterns emerge away from thermodynamic equilibrium, thus thermodynamic theory is not applicable. The understanding of these patterns may require application of kinetic theories.
Pattern dynamics research focuses on universal patternforming mechanisms for the latter case. Away from thermodynamic equilibrium, some experimentally observed 2D spatial patterns have been reported: RayleighBénard convection patterns[8]; reactiondiffusion Turing patterns[9]; Faradaywave patterns[10]; vibratory granular patterns[11]; slimemold spiral patterns[12] and Ca^{2+} spiral waves in Xenopus laevis eggs[13]. Although these patterns are different in the spatiotemporal scales or detailed patternforming 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 reactiondiffusion 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 nonuniform state, leading to the unique textures seen on animal skin. The key patternforming idea in Turing’s paper is the system’s spontaneous symmetrybreaking 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 reactiondiffusion 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 symmetrybreaking 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 patternforming systems essentially consist of coupled differential equations, the simulated patterns are time and spacedependent solutions of the differential equations. The MATLAB builtin fourthorder RungeKutta 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 reactiondiffusion systems. With some modifications, Yang’s programs can be used to simulate pattern formation in a wide range of applications of nonlinear reactiondiffusion 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 scriptbased pattern simulations require highlevel programming skills to tune the model parameters in order to examine the pattern dynamics; this limits their uptake by nonprogrammers. A few attempts for graphicbased pattern simulations have been made in the last decades. For example, READY[33] is an ideal program for getting started with reactiondiffusion systems. It runs fast (taking advantage of multicore architectures), is easy to use (graphicbased) and runs on multiple platforms. In addition, there are other Java applets allowing easy pattern simulations such as GrayScott Simulator[34], showing the GrayScott 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 3D of pattern simulations. However, these graphical user interface (GUI)based softwares designed mainly for demonstrations of specific (preloaded) models. READY supports custom models but is written in READYspecific codes. So technically, READY is still a codingbased 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 biologicallyrelevant patternforming systems.
The aim of this paper is to introduce SIMULINK modelling for simulating pattern dynamics. SIMULINK, an addon 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 modelbased 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 wellknown Brusselator model in SIMULINK. By extending the Brusselator modelling strategies, we construct a SIMULINKbased human cortical model[44] (developed by the authors’ team) that incorporates the patternforming essentials while retaining neurophysiological accuracy: the cortical model comprises interacting populations of excitatory and inhibitory neurons which communicate via chemical (neurotransmittercontrolled) and electrical (gapjunction) 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], deepsleep slowwave oscillations[45] and cognitive gammawaves[46, 47].
Finally, we examine the Brusselator and cortical model pattern dynamics. These simulation results are compared with the theoretical predictions.
Methods
The diffusion constant D _{ U,V } [with units (length)^{2}/time] is an important parameter indicative of the diffusion mobility. For a multicomponent 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 patternforming system in the form of Eq. (1) requires interpreting the differential operators for time ∂/∂ t and space ∇^{2}. 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
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.
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 product block (labelled with ×) in Figure2 combines (1  x ^{2}) 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 realtime output of the two integrators; the scope can be placed anywhere to monitor the response of a subsystem. 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.
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 patternforming 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 differentialequation models in this paper.
Brusselator model
The Brusselator model describes the competition of two chemical species in a chemical reaction, and is the simplest reactiondiffusion 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
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.
In SIMULINK, we initialise the Brusselator model as a column vector consisting of a 60 × 1 grid (spatial resolution = 1 cm/gridpoint) for the onedimensional case; or as a 60 × 60 grid for the twodimensional case. Grid edges are joined to give toroidal boundaries.
where we have again assumed a square grid so that h _{ x } = h _{ y } = h.
The reason we introduce the custom block is that the SIMULINK builtin 2D CONV block provides only the “valid” (nonflux) boundary condition, and cannot handle periodic boundaries.
SIMULINK versions of Waikato cortical model equations
The authors’ team at the University of Waikato (New Zealand) has developed a meanfield 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 (gapjunction) 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.
 1.
Cortical element is the “macrocolumn” containing ∼100,000 neurons.
 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.
There are three isotropic neuronal interactions: corticocortical (longrange from other macrocolumns), intracortical (shortrange within macrocolumn) and connections from subcortical structures such as the thalamus and brainstem. 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 corticocortical connections. Both excitatory and inhibitory connections are permitted in local neuron connections.
 4.
Neurons exchange information via both chemical and electrical (gapjunction) synapses. Gapjunctions are more abundant within inhibitory populations, while being rare within excitatory neuronal populations.
The authors’ team first introduced gapjunctions into the cortical model in the paper Gapjunction mediate largescale Turing structures in a meanfield cortex driven by subcortical noise[44]. The cortical model is expected to exhibit similar dynamics to a chemical reactiondiffusion system: The interaction between the excitatory and inhibitory neurons is analogous to to the competition between the activator and inhibitor of the chemical reactiondiffusion model, e.g. the Brusselator. The gapjunction 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 reactiondiffusion system.
Model equations
where ξ _{ eb } is the Gaussiandistributed whitenoise sources being deltacorrelated in time and space. Inclusion of whitenoise stimulation is motivated by physiological evidence that the cortex requires a continuous background “wash” of input noise to function normally.
We write D _{ bb } as the diffusivecoupling strength between electrically adjoined e ↔ e, i ↔ i neuron pairs. To simplify the notation, we write (D _{1},D _{2}) ≡ (D _{ ee },D _{ ii }).
Symbol definitions for the cortical model
Symbols  Description  Value  Unit 

Λ_{ e,b }  Inverselength 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}$  Longrange e → b axonal connectivity  2000   
${N}_{\mathit{\text{eb}},\mathit{\text{ib}}}^{\beta}$  Local e → b, i → b axonal connectivity  800, 600   
γ _{ e }  Excitatory rateconstant  170  s^{1} 
γ _{ i }  Inhibitory rateconstant  50  s^{1} 
$\u3008{\varphi}_{\mathit{\text{eb}}}^{\text{sc}}\u3009$  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 gapjunction diffusive coupling strength  0–2.0  cm^{2} 
D _{1}  e ↔ e gapjunction diffusive coupling strength  D _{2}/100  cm^{2} 
SIMULINK versions of Waikato cortical model equations
Let us first list the mathematical equations for the Waikato cortical model and examine their characteristics.

The corticocortical 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}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}{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 intracortical 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}\\ \phantom{\rule{.45em}{0ex}}{\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 γ ^{2}Φ to the RHS of the intracortical equations, then the LHS of the intracortical 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 corticocortical 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 rearranged as$\begin{array}{l}\frac{\partial {V}_{b}}{\partial t}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}\frac{1}{{\tau}_{b}}\left[\phantom{\rule{0.3em}{0ex}}{V}_{b}^{\text{rest}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{V}_{b}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}({\rho}_{e}{\psi}_{\mathit{\text{eb}}}{\mathrm{\Phi}}_{\mathit{\text{eb}}}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}{\rho}_{i}{\psi}_{\mathit{\text{ib}}}{\mathrm{\Phi}}_{\mathit{\text{ib}}})\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}{D}_{\mathit{\text{bb}}}{\nabla}^{2}{V}_{b}\phantom{\rule{0.3em}{0ex}}\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).
In following sections, we detail the SIMULINK implementation of the three subsystems (drawn as three blocks in Figure8) of the cortical model.
Corticocortical flux
Intracortical flux
Soma voltage
Results and discussion
Brusselator pattern simulations in 1D space
The onsets of four main classes of instability
Pattern stability  

Bifurcation  Critical wavenumber  Eigenvalue (σ _{ q } = α + iω)  Spatially  Temporally 
Turing  q ≠ 0  σ _{ q } = 0  Unstable  Stable 
Hopf  q = 0  α = 0, ω ≠ 0  Stable  Unstable 
Turing–Hopf  q ≥ 0  α = 0, ω ≠ 0 at q = 0  Unstable  Unstable 
Wave  q ≠ 0  α ≥ 0, ω > 0  Wave instability 
Brusselator pattern simulations in 2D space
Brusselator simulation in a 2D space has more varieties than in a 1D space. These 2D patterns will have unique spatial Turing structures that we cannot infer from the LSA. To precisely predict the Turing pattern dynamics, we derived the amplitude equations[54] for the Brusselator model at the onset of a Turing instability. The amplitude equation describes a reduced form of a reactiondiffusion system yet still retains its essential dynamical features. By approximating the analytic solution, the amplitude equation allows precise predictions of the pattern dynamics when the system is near a bifurcation point. In simple words, the amplitude equations may offer a guidance for model parameter settings corresponding to specific patterns, e.g., Turing pattern textures.
For the 2D spatiotemporal simulation of the Brusselator model, SIMULINK was set to 50s simulation time in auto timestep mode, allowing the patternforming system to organise itself sufficiently. Figure18 predicts the stabilities of stripes (red), H_{0} (blue) and H_{ π } (black) modes for the Turing instability of the Brusselator model. From the range of solid curves, we have the summary of parametric space where a specific mode is stable: The stripes mode is stable when μ s < μ < μ s+; the hexagonal mode is stable (H_{ π } and H_{0}) when μ < μ h or μ > μ h+. H_{ π } and H_{0} interact at μ _{p} where they exchange mode stability, that is, H_{ π } will transit to H_{0} when μ crosses μ _{p} from its LHS to the RHS.
 (a)
μ = 0.0495 falling into a range where only H_{ π } is stable;
 (b)
μ = 0.1100 falling into a range where stripes and H_{ π } modes coexist;
 (c)
μ = 0.3994 where only the stripes structure is stable;
 (d)
μ = 0.7000 again falling into a bistable range where stripes and H_{0} coexist;
 (e)
μ = 1.4802 where only H_{0} is stable.
In Figure18, we see good agreement between SIMULINK simulated patterns and theoretical predictions. Clear H_{ π }, stripes and H_{0} structures are observed at (a), (c) and (e) cases respectively. (b) and (d) show mixed states between forward and backward stable structures. Consequently, we can conclude that by increasing μ, the distance to the critical point, positively, Brusselator forms sequentially from H_{ π }, to stripes, to H_{0}.
In summary, to construct the Brusselator model in SIMULINK, we first place an integrator block to represent time derivative. The temporal integrator’s output will be fed back into the system to engage with the system’s evolution, then form the input of this integrator block, closing the loop. The model parameters can be adjusted by tuning the settings of the blocks A and B as well as two gains (labelled D _{ X } and D _{ Y } respectively). The realtime spatiotemporal evolution of X and Y are monitored via the Matrix viewer block. The simout block delivers the solution of Eq. (6) to the MATLAB workspace for future analysis. The solution is a threedimensional matrix with the third dimension the same length as the time span.
In the next section, we will implement the SIMULINKbased cortical model and examine its clinicallyrelevant pattern dynamics.
Cortical model stability and simulations
At the vicinity of a Turing instability, a weak Hopf instability can be induced in parallel by prolonging the timing of delivery of inhibition at chemical synapses. This permits slow Hopf oscillations with the spatial structure maintained. Specifically, a reduction of the inhibitory rateconstant γ _{ i } in Eq. (20b) below a critical value ∼30.94 s^{1} is sufficient to produce a complex dominant eigenvalue at zero wavenumber whose real part is positive; thus suggesting a global Hopf oscillation.
SteynRoss et al.[46] posit that the interacting lowfrequency Hopf and Turing instabilities may form the substrate for the cognitive state, namely, the “default” background state for the noncognitive brain. These slow patterned oscillations may relate to very slow (≤0.1 Hz) fluctuations in BOLD (bloodoxygenlevel dependent) signals detected using fMRI (functional magnetic resonance imaging) of relaxed, nontasked human brains[62, 63]. SteynRoss et al.[47] also predict that this default state will be suppressed with elevated levels of subcortical drive during goaldirected tasks[64–68].
The effect of embedding MATLAB functions in SIMULINK running efficiency
Although SIMULINK has an intuitive programming logic and comparable accuracy to MATLAB, it sometimes runs much slower than MATLAB, e.g., in the demonstrated Brusselator and cortical model simulations. The reason is that we embed MATLAB functions wraparound in the model to expand the SIMULINK capability. Once a MATLAB function block is present, the MATLAB interpreter is called at each timestep. This drastically reduces the simulation speed. So, one should use the builtin blocks whenever possible. Without using MATLAB function blocks, SIMULINK shows a higher performance than MATLAB, e.g., see the description of Figure3. MathWorks Support Team also presented comprehensive guidance to speed up the SIMULINK simulation, which are available athttp://www.mathworks.com/matlabcentral/answers/94052. In the further optimisation of our SIMULINK model, we will consider replacing MATLAB function with the MEX Sfunction, which may help to accelerate the simulation in the merit of its direct communication with the SIMULINK engine (avoid the time consuming compilelinkexecute cycle). The build of MEX Sfunction requires higher level programming skills, so we only address a more accessible method (embedding MATLAB function) for nonexperts in this paper.
Conclusions
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 codescript and SIMULINK blockdiagram, showing the SIMULINK model has comparable performance with the codescript version.
Using a wellknown Brusselator model, we demonstrated two main SIMULINK modelling strategies for a reactiondiffusion system: interpretation of 2D 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 patternforming 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 SteynRoss et al.’s recent publications[6, 44–47] for comprehensive investigations of the cortical model.
Endnotes
^{a} vanderpoldemo is a MATLAB precoded function
^{b}fixed timestep 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
^{c}The twodimensional 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/22619fast2dconvolution
^{d} MATLAB simulation codes were written by Alistair SteynRoss. The complete codes, plus README files and movies of cortical dynamics, are available from the web site:http://www2.phys.waikato.ac.nz/asr/
Declarations
Acknowledgements
We acknowledge supports from the New Zealand–Japan Exchange Programme (NZJEP), and Embedded Systems Lab in Gunma University, Japan.
Authors’ Affiliations
References
 Ball P: Nature’s Patterns: A Tapestry in Three Parts. 2009, USA: Oxford Univerity PressGoogle Scholar
 Horváth J, Szalai I, De Kepper P: An experimental design method leading to chemical Turing patterns. Science. 2009, 324 (5928): 772775. 10.1126/science.1169973.View ArticlePubMedGoogle Scholar
 Carmon TT, Buljan HH, Segev MM: Spontaneous pattern formation in a cavity with incoherent light. Opt Express. 2004, 12 (15): 34813487. 10.1364/OPEX.12.003481.View ArticlePubMedGoogle Scholar
 Bose S, Rodin P, Scholl E: Competing spatial and temporal instabilities in a globally coupled bistable semiconductor system near a codimensiontwo bifurcation. Phys Rev E. 2000, 62 (2): 17781789. 10.1103/PhysRevE.62.1778.View ArticleGoogle Scholar
 GuiQuan S, Zhen J, QuanXing L, Li L: Pattern formation induced by crossdiffusion in a predator–prey system. Chin Phys B. 2008, 17: 39363941. 10.1088/16741056/17/11/003.View ArticleGoogle Scholar
 SteynRoss ML, SteynRoss DA, Sleigh JW: Gap junctions modulate seizures in a meanfield model of general anesthesia for the cortex. Cogn Neurodynamics. 2012, 6 (3): 215225.View ArticleGoogle Scholar
 Cross MC, Hohenberg PC: Patternformation outside of equilibrium. Rev Mod Phys. 1993, 65 (3): 8511112. 10.1103/RevModPhys.65.851.View ArticleGoogle Scholar
 Gunaratne G, Ouyang Q, Swinney H: Pattern formation in the presence of symmetries. Phys Rev E. 1994, 50 (4): 28022820. 10.1103/PhysRevE.50.2802.View ArticleGoogle Scholar
 Yang L, Dolnik M, Zhabotinsky AM, Epstein IR: Turing patterns beyond hexagons and stripes. Chaos. 2006, 16 (3): 037114037114. 10.1063/1.2214167.View ArticlePubMedGoogle Scholar
 Chen PP, Viñals JJ: Amplitude equation and pattern selection in Faraday waves. Phys Rev E. 1999, 60 (1): 559570.View ArticleGoogle Scholar
 Gollub J, Langer J: Pattern formation in nonequilibrium physics. Rev Mod Phys. 1999, 71 (2): 396403. 10.1103/RevModPhys.71.S396.View ArticleGoogle Scholar
 Nicolis G, Prigogine I: SelfOrganization in Nonequilibrium Systems: From Dissipative Structures to Order Through Fluctuations. 1977, New York: WileyGoogle Scholar
 Wagner J, Li YX, Pearson J, Keizer J: Simulation of the fertilization Ca^{ 2+ } wave in Xenopus laevis eggs. Biophys J. 1998, 75 (4): 208810.1016/S00063495(98)776519.PubMed CentralView ArticlePubMedGoogle Scholar
 Turing AM: The chemical basis of morphogenesis. Philos Trans R Soc Lond B. 1952, 237 (6): 3772.View ArticleGoogle Scholar
 de Wit A: Spatial patterns and spatiotemporal dynamics in chemical systems. Adv Chem Phys. 1999, 109: 435514.Google Scholar
 Garfinkel A, Tintut Y, Petrasek D, Bostrom K, Demer LL: Pattern formation by vascular mesenchymal cells. Proc Natl Acad Sci U S A. 2004, 101 (25): 92479250. 10.1073/pnas.0308436101.PubMed CentralView ArticlePubMedGoogle Scholar
 Rovinsky A, Menzinger M: Selforganization induced by the differential flow of activator and inhibitor. Phys Rev Lett. 1993, 70 (6): 778781. 10.1103/PhysRevLett.70.778.View ArticlePubMedGoogle Scholar
 Jensen O, Pannbacker VO, Mosekilde E, Dewel G, Borckmans P: Localized structures and front propagation in the LengyelEpstein model. Phys Rev E. 1994, 50 (2): 736749. 10.1103/PhysRevE.50.736.View ArticleGoogle Scholar
 Coullet P, Riera C, Tresser C: Stable static localized structures in one dimension. Phys Rev Lett. 2000, 84 (14): 30693072. 10.1103/PhysRevLett.84.3069.View ArticlePubMedGoogle Scholar
 Berezovsky F, Karev G, Song B, CastilloChavez C: A simple epidemic model with surprising dynamics. Math Biosci Eng. 2005, 2 (1): 133152.PubMedGoogle Scholar
 Tlidi M, Mandel P, Lefever R: Localized structures and localized patterns in optical bistability. Phys Rev Lett. 1994, 73 (5): 640643. 10.1103/PhysRevLett.73.640.View ArticlePubMedGoogle Scholar
 Kessler MA, Werner BT: Selforganization of sorted patterned ground. Science. 2003, 299 (5605): 380383. 10.1126/science.1077309.View ArticlePubMedGoogle Scholar
 Maini PK, Woolley TE, Baker RE, Gaffney EA, Lee SS: Turing’s model for biological pattern formation and the robustness problem. Interface Focus. 2012, 2 (4): 487496. 10.1098/rsfs.2011.0113.PubMed CentralView ArticlePubMedGoogle Scholar
 Yang XS: Pattern formation in enzyme inhibition and cooperativity with parallel cellular automata. Parallel Comput. 2004, 30 (5): 741751.View ArticleGoogle Scholar
 Yang XS: An Introduction to Computational Engineering with MATLAB. 2006, London: Cambridge Int Science PublishingGoogle Scholar
 Li J, Chen YT: Computational Partial Differential Equations Using MATLAB. 2011, USA: CRC PressGoogle Scholar
 Garvie MR: Finitedifference schemes for reactiondiffusion equations modeling predatorprey interactions in MATLAB. Bull Math Biol. 2007, 69 (3): 931956. 10.1007/s1153800690623.View ArticlePubMedGoogle Scholar
 Opalska K: Efficient MATLAB simulation of the brusselator. Proc. SPIE, 8903, Photonics Applications in Astronomy, Communications, Industry, and HighEnergy Physics Experiments 2013, 89031U. October 25, 2013, doi:10.1117/12.2035279; [http://dx.doi.org/10.1117/12.2035279],Google Scholar
 Tolle DP, Le Novère N: Meredys, a multicompartment reactiondiffusion simulator using multistate realistic molecular complexes. BMC Syst Biol. 2009, 4: 2424.View ArticleGoogle Scholar
 Goddard NH, Hucka M, Howell F, Cornelis H, Shankar K, Beeman D: Towards NeuroML: model description methods for collaborative modelling in neuroscience. Philos Trans R Soc Lond B Biol Sci. 2001, 356 (1412): 12091228. 10.1098/rstb.2001.0910.PubMed CentralView ArticlePubMedGoogle Scholar
 Vollmer J, Menshykau D, Iber D: Simulating organogenesis in COMSOL: cellbased signaling models. Proceedings of COMSOL Conference 2013. 2013, Rotterdam, 64796479.Google Scholar
 Clauß C, Leitner T, Schneider A, Schwarz P: Objectoriented modelling of physical systems with Modelica using design patterns. System Design Automation. Edited by: Merker R, Schwarz W. 2001, Boston: Springer, 195208.View ArticleGoogle Scholar
 Tim H, Robert M, Andrew T, Tom R, Wills D: Ready. [https://code.google.com/p/reactiondiffusion/],
 Abelson, Adams, Coore, Hanson, Nagpal, Sussman: GrayScott model of reaction diffusion. [http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/],
 Lee KJ, McCormick WD, Ouyang Q, Swinney HL: Pattern formation by interacting chemical fronts. Science. 1993, 261 (5118): 192194. 10.1126/science.261.5118.192.View ArticlePubMedGoogle Scholar
 Lidbeck J: 2D GrayScott reactiondiffusion system. [http://www.aliensaint.com/uo/java/rd/],
 Tyagi AK: MATLAB and SIMULINK for Engineers. 2012, USA: Oxford University PressGoogle Scholar
 Jos van Schijndel AWM: A review of the application of SimuLink Sfunctions to multi domain modelling and building simulation. J Build Perform Simul. 2014, 7 (3): 165178. 10.1080/19401493.2013.804122.View ArticleGoogle Scholar
 Reedy J, Lunzman S: Modelbased design accelerates the development of mechanical locomotive controls. SAE 2010 Commercial Vehicle Engineering Congress. 2010, Chicago: SAE International, 201001–1999.Google Scholar
 Taylor CE, Miller GE: Implementation of an automated peripheral resistance device in a mock circulatory loop with characterization of performance values using Simulink Simscape and parameter estimation. J Med Dev Trans ASME. 2012, 6 (4): 04500110.1115/1.4007458.View ArticleGoogle Scholar
 Graf C, Vath A, Nicoloso N: Modeling of the heat transfer in a portable PEFC system within MATLABSimulink. J Power Sources. 2005, 155 (1): 5259.View ArticleGoogle Scholar
 Chen D, Peng H: A thermodynamic model of membrane humidifiers for PEM fuel cell humidification control. J Dyn Syst Meas Contr. 2005, 127 (3): 424432. 10.1115/1.1978910.View ArticleGoogle Scholar
 Löwe A Agar, D: Chemical Reaction Technology. 2005, Weinheim: WileyVCH Verlag GmbHGoogle Scholar
 SteynRoss ML, SteynRoss DA, Wilson MT, Sleigh JW: Gap junctions mediate largescale Turing structures in a meanfield cortex driven by subcortical noise. Phys Rev E. 2007, 76: 011916011916.View ArticleGoogle Scholar
 SteynRoss ML, SteynRoss DA, Sleigh JW: Interacting TuringHopf instabilities drive symmetrybreaking transitions in a meanfield model of the cortex: a mechanism for the slow oscillation. Phys Rev X. 2013, 3 (2): 021005Google Scholar
 SteynRoss ML, SteynRoss DA, Sleigh JW, Wilson MT: A mechanism for ultraslow oscillations in the cortical default network. Bull Math Biol. 2011, 73 (2): 398416. 10.1007/s1153801095659.View ArticlePubMedGoogle Scholar
 SteynRoss ML, SteynRoss DA, Wilson MT, Sleigh JW: Modeling brain activation patterns for the default and cognitive states. NeuroImage. 2009, 45 (2): 298311. 10.1016/j.neuroimage.2008.11.036.View ArticlePubMedGoogle Scholar
 Cartwright ML: Balthazar van der Pol. J Lond Math Soc. 1960, 35: 367376.View ArticleGoogle Scholar
 Kirschfeld K: The physical basis of alpha waves in the electroencephalogram and the origin of the “Berger effect”. Biol Cybern. 2005, 92 (3): 177185. 10.1007/s0042200505471.View ArticlePubMedGoogle Scholar
 Yamapi R, Filatrella G, AzizAlaoui MA: Global stability analysis of birhythmicity in a selfsustained oscillator. Chaos. 2010, 20 (1): 013114013114. 10.1063/1.3309014.View ArticlePubMedGoogle Scholar
 Peng R, Wang MX: Pattern formation in the Brusselator system. J Math Anal Appl. 2005, 309 (1): 1616.View ArticleGoogle Scholar
 Yu P, Gumel A: Bifurcation and stability analyses for a coupled Brusselator model. J Sound Vib. 2001, 244 (5): 795820. 10.1006/jsvi.2000.3535.View ArticleGoogle Scholar
 Pena B, PerezGarcia C: Stability of Turing patterns in the Brusselator model. Phys Rev E. 2001, 64: 056213056213.View ArticleGoogle Scholar
 Pena B, PerezGarcia C: Selection and competition of Turing patterns. Europhys Lett. 2007, 51 (3): 300306.View ArticleGoogle Scholar
 de Wit A, Lima D, Dewel G, Borckmans P: Spatiotemporal dynamics near a codimensiontwo point. Phys Rev E. 1996, 54: 261View ArticleGoogle Scholar
 Liley DT, Cadusch PJ, Wright JJ: A continuum theory of electrocortical activity. Neurocomputing. 1999, 26: 795800.View ArticleGoogle Scholar
 Rennie CJ, Wright JJ, Robinson PA: Mechanisms of cortical electrical activity and emergence of gamma rhythm. J Theor Biol. 2000, 205 (1): 1735. 10.1006/jtbi.2000.2040.View ArticlePubMedGoogle Scholar
 Robinson PA, Rennie CJ, Wright JJ: Propagation and stability of waves of electrical activity in the cerebral cortex. Phys Rev E. 1997, 56 (1): 826840. 10.1103/PhysRevE.56.826.View ArticleGoogle Scholar
 Cross M, Greenside H: Pattern Formation and Dynamics in Nonequilibrium Systems. 2009, London: Cambridge University PressView ArticleGoogle Scholar
 Kidachi H: On mode interactions in reaction diffusion equation with nearly degenerate bifurcations. Prog Theor Phys. 1980, 63 (4): 11521169. 10.1143/PTP.63.1152.View ArticleGoogle Scholar
 Verdasca J, de Wit A, Dewel G, Borckmans P: Reentrant hexagonal Turing structures. Phys Lett A. 1992, 168 (3): 194198. 10.1016/03759601(92)905746.View ArticleGoogle Scholar
 Fox MDM, Snyder AZA, Vincent JLJ, Corbetta MM, Van Essen DCD, Raichle MEM: The human brain is intrinsically organized into dynamic, anticorrelated functional networks. PNAS. 2005, 102 (27): 96739678. 10.1073/pnas.0504136102.PubMed CentralView ArticlePubMedGoogle Scholar
 Fransson P: Spontaneous lowfrequency BOLD signal fluctuations: an fMRI investigation of the restingstate default mode of brain function hypothesis. Hum Brain Mapp. 2005, 26 (1): 1529. 10.1002/hbm.20113.View ArticlePubMedGoogle Scholar
 Raichle ME, MacLeod AM, Snyder AZ, Powers WJ, Gusnard DA, Shulman GL: A default mode of brain function. PNAS. 2001, 98 (2): 676682. 10.1073/pnas.98.2.676.PubMed CentralView ArticlePubMedGoogle Scholar
 Laufs HH, Kleinschmidt AA, Beyerle AA, Eger EE, SalekHaddadi AA, Preibisch CC, Krakow KK: EEGcorrelated fMRI of human alpha activity. NeuroImage. 2003, 19 (4): 1414.View ArticleGoogle Scholar
 Gusnard DA, Raichle ME: Searching for a baseline: Functional imaging and the resting human brain. Nat Rev Neurosci. 2001, 2 (10): 685694. 10.1038/35094500.View ArticlePubMedGoogle Scholar
 Greicius MDM, Krasnow BB, Reiss ALA, Menon VV: Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. PNAS. 2003, 100 (1): 253258. 10.1073/pnas.0135058100.PubMed CentralView ArticlePubMedGoogle Scholar
 Damoiseaux JS, Rombouts SARB, Barkhof F, Scheltens P, Stam CJ, Smith SM, Beckmann CF: Consistent restingstate networks across healthy subjects. PNAS. 2006, 103 (37): 1384813853. 10.1073/pnas.0601417103.PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.