# Stable moving patterns in the 1-D and 2-D Gray-Scott Reaction-Diffusion System

I presented some of the material from my 2009 paper, and a few of the more interesting videos from the xmorphia exhibit to the Mathematics department at Rutgers University as part of their Fall 2010 series in Mathematical Physics.

Here are the slides in two formats, Powerpoint and Adobe PDF:

Stable moving patterns (Rutgers talk) -- slides (PPS) (13.1 MB)

Stable moving patterns (Rutgers talk) -- slides (PDF) (12.6 MB)

The slides in HTML format are below (scroll down past the five large color images)

The PPS and PDF versions contain extensive graphics and links to the videos I prepared for the talk (which are on YouTube). Some software (e.g. PDF plug-in or a PPS viewer) might not be able to jump to the links directly, but you can cut and paste the URL text (example: "youtube.com/v/kXDTqqgrYCg") into your web browser's address bar.(The same links are also in the slide text below)

In addition to the videos linked from the slides, the participants also saw the videos from the following pages in my main Gray-Scott web exhibit:

*k=0.055, F=0.098*

(showed briefly while looking for a better example)

*k=0.055, F=0.102*

(showed briefly while looking for a better example)

*k=0.057, F=0.098*

(persistent stripes that form a linked network of polygonal cells,
with gradual evolution into fewer and larger cells reminiscent of soap
bubbles)

*k=0.057, F=0.014*

(showed briefly while looking for a better example)

*k=0.055, F=0.014*

(low-U high-V spots that move a bit, then "split" into two spots
which then move and split again, etc.)

What follows is the text of each slide with links to the relevant images and videos. This is intended for readers with a slower connection, or who prefer to go to the YouTube videos directly from the browser without the need to download the slides.

Stable moving patterns in the 1-D and 2-D

Gray-Scott Reaction-Diffusion System

*figure: two spots far*

*figure: U-shaped pattern*

*figure: 4 spots zigzag*

*figure: three spots 1-D*

*figure: five spots*

*figure: "landspeeder"*

Robert Munafo

Rutgers Mathematical Physics Seminar Series

2010 Dec 9, 2 PM — Hill 705

slide 1

Contents (outline of this talk)

- Brief history, motivation and discovery

- Methods and testing

- New patterns and interactions (mostly video)

- Connections; Open questions

- Discussion

slide 2

Brief History: Initial Motivation

- 1994: Looking for a problem to run on new hardware
- Supercomputer research led to Williams exhibit at Caltech (shown, left)
- Literature was easy to find; problem was appealing
- Exploring parameter space more closely

*figure: 1994 Williams Caltech website*

Roy Williams "Xmorphia" web exhibit (Caltech, 1994)

*figure: PDE formulas*

Pearson "Complex patterns in a simple system" (Science 261 1993) (illust. next slide)

Lee et al. "Experimental observation of self-replicating spots in a reaction-diffusion system" (Nature 369 1994)

(showing part of fig. 1 from Lee et al., and their form of the PDE formulas)

Ferrocyanide-iodate-sulphite reaction in gel reactor

( K_{4}Fe(CN)_{6}^{.}3H_{2}O , NaIO_{3} , Na_{2}SO_{3} , H_{2}SO_{4} ,
NaOH , bromothymol blue )

slide 3

Brief Intro: Parameter Space

Parameter space and four pattern examples, from Pearson (ibid. 1993)

(showing fig. 3 and part of fig. 2 from Pearson (ibid. 1993))

Coloring imitates bromothymol blue pH indicator as in Lee et al. "Pattern formation by interacting chemical fronts" (Science 261 1993): blue = low U; yellow = intermediate; red = high U

(showing fig. 3 from Miyazawa et al.)

Parameter space visualizations for several reaction- diffusion models, from fig. 3, Miyazawa et al. "Blending of animal colour patterns by hybridization" (Nature Comm. 1(6) 2010)

slide 4

Brief History: 2009 Website Project

*figure 1: parameter map image from my website*

- Denser coverage; explore more extreme parameter values

Pearson: 34 sites; 0.045

- Several types of starting patterns

- Video for each (k,F) point

- Higher resolution and precision

- Coloring to show both U and dU/dt

- Describe and catalog all phenomena

- Run each pattern "to completion" no matter how long that takes

slide 5

New pattern types (2009)

*negative soap-bubbles image*

k=0.057, F=0.09

*positive soap-bubbles image*

k=0.059, F=0.09

*exotic diversity image*

k=0.061, F=0.062

*self-sustaining spirals image*

k=0.045, F=0.014

(The above 4 images have arrows pointing to a copy of the parameter map image from the previous slide, and showing the locations of the k and F values for each image)

- Some type of exotic behavior at high F values was expected

- Spirals were expected (seen in Belousov-Zhabotinsky and other reaction-diffusion systems)

- k=0.061, F=0.062 has more diversity than all the rest put together

- Also found: mixed spots and stripes; variations in branching; etc.

(Xmorphia parameter map image again)

Munafo (2009) www.mrob.com/pub/comp/xmorphia

slide 6

Inherent Instability

Turing-F600-k620.mp4 — watch (YouTube)
— youtube.com/v/kXDTqqgrYCg

k=0.062, F=0.06

Periodic boundary conditions; size 3.35w × 2.33h

Each second is 1102 dtu

Initial pattern of low-level random noise (0.4559<U<0.4562; 0.2674<V<0.2676)

Final values: 0.35<U<0.90; 0.00<V<0.36 (approx.)

(Contrast-enhanced images; lighter = higher U)

- At many parameter values, patterns like this grow out of any inhomogeneity, no matter how small

- This is not a consequence of numerical approximation error: proven by mathematical analysis. Dominant wavelength (spot size) depends on reaction dynamics and diffusion rate. Turing "The chemical basis of morphogenesis" (Phil. Trans. Royal Soc. London B 237(641) 1952)

*figure: parameter space with Turing and exotic regions*

Gray-Scott parameter values for Turing instability (green; width to
scale) and stable moving patterns (blue; width exaggerated). The video
is at X

slide 7

How can we trust these results?

- Arbitrarily low-level noise can generate strong patterns

- Instabilities can exist in the ideal (exact) system, and can be introduced by the numerical method

- Complex patterns are already proven to be real

- These new patterns are far more extraordinary

- Mathematical proof/disproof probably impossible

- How much precision is necessary? Is any finite precision sufficient for proof?

- Is the standard precision "too accurate" to be relevant? The real world has known levels of quantization and randomness: "finite precision"

- Goals:

- Eliminate suspicion of numerical error

- Quantify the sensitivity of these pattern phenomena to precision, randomness, reaction parameters and other environmental conditions

slide 8

Supplement to Slide 8

My diagram illustrating the effects of "precision", "noise" and "randomness" on the robustness of pattern-forming phenomena

^ Q | R L | | ^ | : | : | : | : | M MO:M OM P ---E--------:-----------> | NThis is an abstract (qualitative conceptual) representation of two different types of phenomena that affect pattern-forming systems, viewed as "dimensions" or "parameters" of the laws of physics.

Q represents the amount of quantization that results from the (finite, not infinitesimal) size of atoms. This determines the relative strength of brownian motion and other similar phenomena.

N represents the amount of aberration resulting from numerical methods — roundoff error, the choice of a certain finite grid size to represent a continuous system, etc.

The other letters represent the location of certain pattern-forming systems that have been studied in published work:

P : Pearson (ibid. 1993)

M : Munafo (my various Gray-Scott simulation experiments from 1994 to 2009)

O : Other published Gray-Scott simulations over the period 1994-2009

E : Exact systems (mathematical analysis and proofs)

R : The real world

L : Stochastic simulations, such as certain papers by Prof. Joel Lebowitz of the Rutgers mathematics department.

The vertical dotted line represents one author (Munteanu et al. "Pattern formation in noisy self-replicating spots" Int'l. Jrl. of Bifurcation and Chaos 12(16) 2006) who added various amounts of "Gaussian noise" to a simulation and observed a substantial qualitative change in pattern type (from stripes to spots) as one goes from zero noise up to a level one might expect in real-world experiments, in Gray-Scott simulations with k=0.0655, F=0.05 and all other details as in Pearson (ibid. 1993). I have reproduced these results in my simulations.

Note that these "systems" exist in different types of "reality" : the physical reality of the world around us, the "virtual reality" inside a computer, and the "abstract" existence of mathematical truth. For this and other reasons, there are many more than just two "dimensions" at play and so the chart above is just an illustration of the concept.

I sought to address the general issue of numerical error ("M" systems as compared to "E") and ignoring statistical mechanics ("E" systems as compared to "R").

Verification Examples

*figure: two spots (close)*

*figure: U-shaped pattern*

- Two stable moving phenomena

- One is clearly bogus, the other might be real

Two spots maintain the same distance while the pair rotates toward
a 45^{o} alignment

U-shaped pattern moves at about 1 dlu per 62,000 dtu (dimensionless units of length, time)

- Define something that is measurable

- Model the sources of error, e.g.:

measured value = true value + simulation error + measurement error

simulation error = f(stability, precision, grid spacing, ...)

CFL (Courant-Friedrichs-Lewy 1928) stability criterion (for the Laplacian term):

C Δx^{2}/Δt

(A higher value means greater stability. Constant C depends on e.g. k and F for the Gray-Scott system)

- Progressively improve the simulation and look for a trend

slide 9

Verification Procedure

*screen shot of the same test running in five different-sized grids*

- Progressively improve Δx by √2

- Progressively improve CFL stability by √2

- Progressively improve Δt by needed amount (2√2)

Δx | Δt | Δx^{2}/Δt | pixels(typ.) | |

"std" | 0.00699 | 0.5 | 9.78×10^{-5} | 128×128 |

"s1.4" | 0.00495 | 0.177 | 1.38×10^{-4} | 180×180 |

"s2" | 0.00350 | 0.0625 | 1.96×10^{-4} | 256×256 |

"s2.8" | 0.00247 | 0.0221 | 2.77×10^{-4} | 360×360 |

"s4" | 0.00175 | 0.00781 | 3.91×10^{-4} | 512×512 |

(amount of calculation increases by 4√2 each time: ratio of 1024 to 1 between "std" and "s4")

Expected Results:

(two-spot pattern): Movement is 100% spurious: measurements should tend towards zero

(U-shaped pattern): If real, velocity should clearly converge on a nonzero value

slide 10

Verification Examples — Results

*figure: two spots (close)*

Suspect rotating 2-spot phenomenon

cross-qtr. relative model max dU/dt interval velocity std 8.97e-7 5.6e5 1.0 s1.4 4.57e-7 1.12e6 0.50 s2 2.31e-7 2.18e6 0.26 s2.8 1.158e-7 4.34e6 0.129 s4 5.80e-8 8.7e6 0.064*figure: U-shaped pattern*

Velocity of U-shaped pattern

model dist pixels time velocity (meas.err) std 0.55418 79.2 35076 1/63294(45) s1.4 0.59863 121.1 37342 1/62379(64) s2 0.57042 163.1 35456 1/62159(27) s2.8 0.56947 230.3 35258 1/61913(44) s4 0.56553 323.5 35009 1/61905(12)Same tests using 4th order Runge-Kutta

model dist pixels time velocity (meas.err) std 0.57937 82.9 36559 1/63101 s1.4 0.56053 113.4 35021 1/62479 same s2 0.56351 161.2 35008 1/62124 as s2.8 0.56517 228.6 35003 1/61934 above s4 0.56574 323.6 35002 1/61870Limits on parameter k (when F=0.06) for stability of U-shaped pattern

model minimum maximum std 0.0608833 0.0609829 s1.4 0.0608796 0.0609831 s2 0.0608777 0.0609831 s2.8 0.0608767 0.0609830 s4 0.0608762 did not test measurement error in all values is +-1 in the last digit- Two-spot pattern movement is bogus

- Moving U-shaped pattern is real; Runge-Kutta gives little if any benefit

- Asymptotic trends in values and in measurement errors, as expected

slide 11

Immunity to Noise

U-noise-immunity.mp4 — watch (YouTube)
— youtube.com/v/_sir7yMLvIo

k=0.0609, F=0.06

Periodic boundary conditions; size 3.35w × 2.33h

1 second ≈ 1100 dtu

4 U-shaped patterns traveling "up"

Systematic noise perturbation applied once per 73.5 dtu

Amplitude of each noise event starts at 0.001 and doubles every 11,000 dtu (10 seconds in this movie)

At noise level 0.064, three patterns are destroyed; noise amplitude is then diminished to initial level

(False-color: yellow = high U; pastel = positive dU/dt )

- Pattern continues to move in the presence of noise, and generally behaves as expected of a real phenomenon

- Similar tests (e.g. more frequent noise events each of lower amplitude) give similar results

slide 12

Symmetry-based Instability Tests

Daedalus-stability.mp4 — watch (YouTube)
— youtube.com/v/fWfsMVEeP5k

k=0.0609, F=0.06

Periodic boundary conditions; size 2.36w × 1.65h

Manually constructed initial pattern based on parts of naturally-evolved systems

1 second = 265 dtu

Two sets of three noise events; noise event amplitudes are 10^{-5},
0.001 and 0.1; pattern allowed to recover after each event

Coloring (left image) same as before

Coloring (images below): white = positive dU/dt, black =
negative dU/dt, shades of gray when magnitude of dU/dt is less
than about 5×10^{-13}

*figure: landspeeder (full color)*

*figure: landspeeder (derivative only)*

A pattern with an instability that this test does not help reveal

- Rotating and moving patterns have rotational or bilateral (resp.) symmetry which persists if the pattern is stable.

- Instability and/or a return to symmetry are easier to observe in the derivative

- Test static and linear-moving patterns at different angles to reveal influence by grid effects

- Other tests include shifting parts of a pattern, applying noise to only part of a pattern, etc.

- Such tests can reveal instability more quickly but do not prove stability.

slide 13

Logarithmic Timebase, Long Duration

- Some phenomena appear asymptotic to stability but actually keep moving forever

- Running a simulation for as long as possible (currently >10
^{8}time units) and viewing the results at an exponentially accelerating speed can reveal some of these phenomena

- Repulsion of solitons in 1-D (pictured, lower-right) has been
studied mathematically for high ratios D
_{U}/D_{V}(Doelman et al. "Slowly modulated two-pulse solutions...", SIAM Jrl. on Appl. Math. 61(3) 2000) with the result: (speed of movement) = Ae^{-Bt}for positive constants A,B

- There are many other phenomena in Gray-Scott systems that slow at similar rates

Eight solitons in a 2-D system, k=0.067, F=0.046 :

*t=78,125*

;

*t=1,250,000*

;

*t=2×10 ^{7}*

;

*t=3.2×10*

^{8}*figure 2: pairs of solitons in 2-D and 1-D*

(larger, upper part of figure) Pair of solitons in 2-D system,
k=0.068, F=0.042

(smaller, lower part of figure) Pair of
solitons in 1-D system, k=0.0615, F=0.04

slide 14

Discovery - Great Diversity

Original-F620-k610-fr159.mp4 — watch (YouTube)
— youtube.com/v/wFtXwFfrwWk

k=0.061, F=0.062

Periodic boundary conditions; size 3.35w × 2.33h

1 second ≈ 1190 dtu

Initial pattern of a few randomly placed spots of relatively high U on a "blue" background (secondary homogeneous state, approx. U=0.420, V=0.293)

(False-color: yellow = high U; pastel = positive dU/dt )

*figure 1: normal soliton*

Ordinary spot soliton (k=0.067, F=0.062)

*figure 2: negative soliton*

"negative soliton" (k=0.061, F=0.062)

In figures 1 and 2: White curve = U; Black curve = V; dotted line shows cross-section taken

- Part of routine scan for website exhibit project

- "Negative solitons" (hereafter called "negatons") exhibit attraction and multi-spot binding

- Several types of patterns in one system

- The "target" pattern is not stable at these parameter values, but is stable at the nearby parameters k=0.0609, F=0.06

("iota" image from figure 2 in Pearson (ibid. 1993))

"negative solitons" in Pearson (ibid. 1993)

"Pattern ι is time independent and was observed for only a single parameter value."

(parameters unpublished, probably k=0.06, F=0.05)

slide 15

Discovery - Moving U Pattern

U-discovery-F620-k609-fr521.mp4 — watch (YouTube)
— youtube.com/v/xGMuuPYhLiQ

Original coloring — youtube.com/v/ypYFUGiR51c

k=0.0609, F=0.062

Periodic boundary conditions; size 3.35w × 2.33h

1 second ≈ 3900 dtu

Initial pattern and coloring same as before

*figure: parameter space with ι and X*

This video is at X; Pearson's type ι is shown

- Moving "U" visible to right of center, short-lived (hits two negatons)

- More unexpected negaton behavior: being "dragged" by other features

slide 16

Long Duration Test

Different Behaviors at Multiple Time Scales

Exponential-time-lapse.mp4 — watch (YouTube)
— youtube.com/v/-k98XOu7pC8

k=0.0609, F=0.062

Periodic boundary conditions; size 3.35w × 2.33h

Initial pattern of spots (randomly chosen U<1, V>0) on solid red (U=1, V=0) background; coloring same as before

Video uses accelerating time-lapse: simulation speed doubles every 6.7 seconds

- First 30,000 dtu (40 seconds into video): blue spots grow to fill the space (a very common Gray-Scott behavior)

- Up to 1.25 million dtu (75 seconds into video): complex behavior, double-spaced stripes, solitary negatons, etc. until all empty space is filled

- Up to 6 million dtu (90 seconds into video): chaotic oscillation of parallel stripes growing, shrinking and twisting; gradually producing more spots

- Sudden onset of stability: All chaotic motion ends (whole system drifts very slowly)

slide 17

Complex Interactions

complex-interactions-1.mp4 — watch (YouTube)
— youtube.com/v/hgTBOf7gg8E

k=0.0609, F=0.06

Periodic boundary conditions; size 3.35w × 2.33h

Each second is 10,000 dtu

Manually constructed initial pattern based on parts of naturally-evolved systems

Coloring same as before

- U-shape can influence other objects and survive (although it generally does not)

- The clusters of negatons along the top move and rotate, very slowly

slide 18

Slow Movement, Rotation

slow movers and rotaters.mp4 — watch (YouTube)
— youtube.com/v/PB3lPMhwIo0

k=0.0609, F=0.06

Periodic boundary conditions; size 3.35w × 2.33h

Each second is 100,000 dtu

Manually constructed initial pattern based on parts of naturally-evolved systems

Coloring same as before

- There are many very slow patterns; a few of the more common are shown here

- The "target" (negaton with annulus) is a common product of spotlike initial patterns, and adjacent negatons typically make it move or rotate

- The 4 negatons left by the collision of two U patterns are in an unstable equilibrium

*figure 1: rotating "target with 9 spots" undergoing tests*

Stability analysis of another slow-rotating pattern

*figure 2: 4-spot pattern undergoing tests*

Stability testing of 4-negaton pattern; stable form shown at right

slide 19

A Gray-Scott Pattern Bestiary

*figure: a large collection of patterns*

k=0.0609, F=0.06 Contrast-enhanced grayscale, lighter = higher U

Most patterns were manually constructed from parts of naturally occurring forms

Note: Many of these are not yet thoroughly tested

- Almost everything that keeps its shape moves indefinitely. In general, a pattern will:

- move in a straight line, if is has (only) bilateral symmetry

- rotate, if it has (only) rotational symmetry

- move on a curved path, if it has no symmetry at all

- Lone negatons are frequently "captured" and/or "pushed" by a moving pattern

slide 20

Relation to Other Work:

"Negaton" Clusters and Targets

- Spots and "target" patterns very similar to the Gray-Scott "negatons" are seen in papers by Schenk, Purwins, et al.

- In a 1998 work, a 2-component R-D system is studied; the spots are stationary and are reported to "bind" into stable groupings with specific geometrical configurations (as shown). There are differences between this system and Gray-Scott, evident in which "molecules" are reported as stable.

- In a 1999 work (by Schenk alone) a 3- component R-D system includes a "target" pattern with a very similar cross- section.

(Showing parts of fig. 1 and fig. 5 from Schenk et al. (1998))

From Schenk et al. "Interaction of self-organized quasiparticles..." (Physical Review E 57(6) 1998), fig. 1 and 5 The pattern marked * (subfigure "b") is stable in Gray-Scott at {k=0.0609, F=0.06}, the pattern marked × (subfigure "e") is not.

*figure: target pattern*

Target pattern in Gray-Scott, k=0.0609, F=0.06, with U and V
levels at cross-section through center

(Showing part of fig. 4.37 from Schenk PhD dissertation)

Target pattern from C. P. Schenk (PhD dissertation, WWU Münster, 1999) p. 116 fig. 4.37

slide 21

Relation to Other Work: Halos

- The light and dark "rings" or "halos" are seen in physical experiments and numerical simulations intended to model both physical and biological systems

- When spots appear in these systems, as in Gray-Scott, the spots tend to be seen at certain "quantized" distances

- Halo amplitudes, spot spacings and relative sizes differ; this also reflects changes seen in Gray-Scott as the k and/or F parameters are changed

*figure: two "negatons" at stable weakly-bound distance*

Negatons with halos (Gray-Scott system, k=0.0609, F=0.06, lighter = higher U; exaggerated contrast)

(showing fig. 11 from Barrio et al. 2009)

Spots with halos in numerical simulation by Barrio et al. "Modeling the skin patterns of fishes" (Physical Review E 79 031908 2009) fig. 11

(showing fig. 3d from Stollenwerk web page)

Spots with halos in gas-discharge experiment by Lars Stollenwerk ("Pattern formation in AC gas discharge systems", website of the Institute of Applied Physics, WWU Münster, 2008) fig. 3d

Halos also appear in models by Schenk, Purwins, et al (ibid., 1998 and 1999, shown on previous slide) in work related to gas discharge experiments

slide 22

One-Dimensional Gray-Scott Model

- There are extensive results on the 1-D system based on rigorous
mathematical analysis (most are for higher ratios D
_{U}/D_{V}than in the systems presented here)

- The "spiral wavefront" observed at many parameter values in the 2-D system is also a viable self- sustaining stable moving pattern in the 1-D system at the same parameter values

- For many parameter values that support "negative stripes" in the 2-D system, certain asymmetrical clumps of "negatons" form stable moving patterns in the 1-D system

- Negatons in 1-D were reported (shown) as early as 1996

(showing fig. 9 from Mazin et al. 1996)

Single 1-D negaton inside a growing region of solid "blue state" at k=0.06, F=0.05, from Mazin et al. "Pattern formation in the bistable Gray-Scott model" (Math. and Comp. in Simulation 40 1996) fig. 9

*figure 1: spiral and fast pulse*

2-D spirals and 1-D pulse at k=0.047, F=0.014

*figure 2: two slow-moving patterns*

Representative 2-D and 1-D patterns at k=0.0609, F=0.06

(Note: Non-existence results of Doelman, Kaper and Zegeling (1997)
and of Muratov and Osipov (2000) are not applicable because they
concern models with a much higher ratio D_{U}/D_{V} )

slide 23

Open Questions

- Why does the U-shaped pattern move and keep its shape?

— As parameter k is increased, leading end of double-stripe (shown) moves faster, but trailing end moves slower and the object lengthens; in the other direction (decreasing k) the reverse is true

— When these two speeds are closely matched, the U shape (shown) neither grows nor shrinks — why?

- Do these patterns appear in other reaction-diffusion models?

— Universal presence of other pattern types suggests this; parameter space maps should make it easy to find; nearby Turing effect is possibly relevant

- Can any of the special properties of these patterns be proven mathematically?

— 1-D systems seem particularly well suited to this task

— Shape of negaton "halos" is easy to solve

— Most existing work applies conditions or limits that exclude the
commonly studied D_{U}/D_{V}=2 systems

slide 24

Discussion

*screen image of my Gray-Scott program in July 1994*

18 MFLOPs/sec. Note ASCII text rendering of pattern

*my Gray-Scott video screensaver, September 2010*

(depicts a Gray-Scott image using a brightness map of my face to
control parameter k)

Robert Munafo mrob.com/sci

(email: click "contact" below)

(the background of this slide is a detailed parameter map image generated in 2004, showing k from 0.0582 to 0.0682 and F from 0.031 to 0.051, in a color scheme that depicts low U values as blue and high U values as yellow.)

slide 25

This page was written in the "embarrassingly readable" markup language RHTF, and was last updated on 2019 Jan 05. s.11