Belousov-Zhabotinsky reaction

Whilst exploring Wikipedia I came across the Belousov-Zhabotinsky reaction, and the phenomenology of its reaction dynamics such as the development of propagating waves of reactants is described here. Although I previously knew about the BZ reaction I had never implemented a simulation of it. It turns out to be remarkably easy to implement in Mathematica, and the BZ reaction dynamics can readily be displayed as a very pretty movie.

Although the original purpose of this type of simulation was to model the BZ reaction dynamics, you could imagine using variants of this type of simulation for generating interesting dynamic art forms.

A slightly modified form of an algorithm to simulate the BZ reaction due to Professor A K Dewdney is listed here.

(i) Select an integer q in the range 2 through 255. Cells may be in any of the states 1 through q.
(ii) Select two integers k1 and k2 in the range 1 through 8 and an integer g in the range 0 through 100.
(iii) In the transition from one "step" to the next the state of each cell is changed once according to rules (iv) - (vii).
(iv) A cell in state q changes to state 1.
(v) A cell in state 1 changes to state a/k1 + b/k2 + 1 where a is the number of neighbors of the cell which are in states 2 through q-1 and b is the number of neighbors in state q.
(vi) A cell in any of states 2 through q-1 changes to S/(9 - c) + g, where S is the sum of the states of the cell and its neighbors and c is the number of neighbors in state 1.
(vii) If the application of rule (v) or rule (vi) would result in a cell having a state > q then the state of that cell becomes q.

After some experimentation to find the fastest (and clearest) implementation in Mathematica, I created a function for implementing a BZ update step with toroidal boundary conditions:

update[state_, q_, k1_, k2_, g_] :=
Module[
{ones, qs, kernel, total, totalones, totalqs},
ones = Map[If[#==1, 1, 0]&, state, {2}];
qs = Map[If[#==q, 1, 0]&, state, {2}];
kernel = Table[1, {3}, {3}];
total = ListCorrelate[kernel, state, {{2,2},{2,2}}];
totalones = ListCorrelate[kernel, ones, {{2,2},{2,2}}] - ones;
totalqs = ListCorrelate[kernel, qs, {{2,2},{2,2}}] - qs;
MapThread[If[#>q, q, #]&[Switch[#1, 1, Floor[#4/k1+#3/k2+1], q, 1, _, Floor[#5/(9-#2)+g]]]&, {state, totalones, totalqs, 8-totalones-totalqs, total}, 2]
];

I used this BZ update function to generate the above movie of the BZ reaction dynamics, using suitable parameter values that I copied from here.

statesequence =
With[
{n=64, q=200, k1=2, k2=3, g=70},
NestList[update[#, q, k1, k2, g]&, RandomInteger[{1,q}, {n,Floor[4/3n]}], 250]
];
movie = Map[ArrayPlot[#, ColorFunction->"Rainbow"]&, statesequence];
ListAnimate[movie]

The poor image quality in the movie is due to the high degree of image compression, rather than due to any limitations of the above algorithm. Also, each frame is separately scaled and mapped onto the colour lookup table, which is probably not the best way of creating the movie.