User Tools

Site Tools


examples:reaction-diffusion

Reaction-diffusion models

1D reaction-diffusion: Activator-Inhibitor

Space-time plot of 1D reaction diffusion model.

Introduction

The first example models a 1D activator-inhibitor model (Meinhardt and Gierer, 1972).

Model description

This 1D PDE model uses a Lattice with linear structure and periodic boundary conditions.

The PDE defined two species called Layers: A (activator) and I (inhibitor) with resp. low and high Diffusion rates. The reaction part of the equations are defined in the System. Similar to ODE model, the System defines Constants and DiffEqns.

The results are recorded and visualized using the Logger and SpaceTimeLogger.

Model

ActivatorInhibitor_1D.xml

<MorpheusModel version="3">
    <Description>
        <Title>Example-ActivatorInhibitor1D</Title>
        <Details></Details>
    </Description>
    <Global>
        <Field symbol="a" value="rand_norm(0.5,0.1)" name="activator">
            <Diffusion rate="0.02"/>
        </Field>
        <Field symbol="i" value="0.1" name="inhibitor">
            <Diffusion rate="1"/>
        </Field>
        <System solver="runge-kutta" time-step="1">
            <DiffEqn symbol-ref="a">
                <Expression>(rho*a^2) / i - mu_a * a + rho_a</Expression>
            </DiffEqn>
            <DiffEqn symbol-ref="i">
                <Expression>(rho*a^2) - mu_i * i</Expression>
            </DiffEqn>
            <Constant symbol="rho_a" value="0.01"/>
            <Constant symbol="mu_i" value="0.03"/>
            <Constant symbol="mu_a" value="0.02"/>
            <Constant symbol="rho" value="0.001"/>
        </System>
    </Global>
    <Space>
        <Lattice class="linear">
            <Size symbol="size" value="100 0 0"/>
            <BoundaryConditions>
                <Condition boundary="x" type="periodic"/>
            </BoundaryConditions>
            <NodeLength value="0.25"/>
            <Neighborhood>
                <Order>1</Order>
            </Neighborhood>
        </Lattice>
        <SpaceSymbol symbol="space"/>
    </Space>
    <Time>
        <StartTime value="0"/>
        <StopTime value="4000"/>
        <SaveInterval value="0"/>
        <!--    <Disabled>
        <RandomSeed value="1"/>
    </Disabled>
-->
        <TimeSymbol symbol="time"/>
    </Time>
    <Analysis>
        <Logger time-step="25">
            <Input>
                <Symbol symbol-ref="a"/>
                <Symbol symbol-ref="i"/>
            </Input>
            <Output>
                <TextOutput/>
            </Output>
            <Plots>
                <Plot time-step="250" name="space plot">
                    <Style style="lines" line-width="3.0"/>
                    <Terminal terminal="png"/>
                    <X-axis>
                        <Symbol symbol-ref="space.x"/>
                    </X-axis>
                    <Y-axis minimum="0" maximum="3.5">
                        <Symbol symbol-ref="a"/>
                        <Symbol symbol-ref="i"/>
                    </Y-axis>
                    <Range>
                        <Time mode="current"/>
                    </Range>
                </Plot>
                <Plot time-step="-1" name="time-space plot">
                    <Style style="points"/>
                    <Terminal terminal="png"/>
                    <X-axis>
                        <Symbol symbol-ref="time"/>
                    </X-axis>
                    <Y-axis>
                        <Symbol symbol-ref="space.x"/>
                    </Y-axis>
                    <Color-bar>
                        <Symbol symbol-ref="a"/>
                    </Color-bar>
                </Plot>
            </Plots>
        </Logger>
    </Analysis>
</MorpheusModel>


In Morpheus GUI: Examples → PDE → ActivatorInhibitor_1D.xml


2D reaction-diffusion: Activator-Inhibitor

Stripe pattern generated by 2D Gierer-Meinhardt model

Introduction

A 2D activator-inhibitor model (Meinhardt and Gierer, 1972).

Description

This model uses a standard Lattice with square structure and periodic boundary conditions.

The PDE defined two species or Layers A (activator) and I (inhibitor) with resp. low and high Diffusion rates. The reaction part of the equations are defined in the System.

The results are visualized using the Gnuplotter.

Model

ActivatorInhibitor_2D.xml

<MorpheusModel version="3">
    <Description>
        <Title>Example-ActivatorInhibitor-2D</Title>
        <Details></Details>
    </Description>
    <Global>
        <Field symbol="a" value="rand_norm(0.5,0.1)" name="activator">
            <Diffusion rate="0.02"/>
        </Field>
        <Field symbol="i" value="0.1" name="inhibitor">
            <Diffusion rate="1"/>
        </Field>
        <System solver="runge-kutta" time-step="5" name="Meinhardt">
            <Constant symbol="rho" value="0.001"/>
            <Constant symbol="rho_a" value="0.001"/>
            <Constant symbol="mu_i" value="0.03"/>
            <Constant symbol="mu_a" value="0.02"/>
            <Constant symbol="kappa" value="0.10"/>
            <DiffEqn symbol-ref="a">
                <Expression>(rho/i)*((a^2)/(1 + kappa*a^2)) - mu_a * a + rho_a</Expression>
            </DiffEqn>
            <DiffEqn symbol-ref="i">
                <Expression>rho*((a^2)/(1+kappa*a^2)) - mu_i *i</Expression>
            </DiffEqn>
        </System>
    </Global>
    <Space>
        <Lattice class="square">
            <Size symbol="size" value="150 150 0"/>
            <NodeLength value="1"/>
            <BoundaryConditions>
                <Condition boundary="x" type="periodic"/>
                <Condition boundary="y" type="periodic"/>
            </BoundaryConditions>
            <Neighborhood>
                <Order>1</Order>
            </Neighborhood>
        </Lattice>
        <SpaceSymbol symbol="space"/>
    </Space>
    <Time>
        <StartTime value="0"/>
        <StopTime value="12500"/>
        <SaveInterval value="0"/>
        <RandomSeed value="2"/>
        <TimeSymbol symbol="time"/>
    </Time>
    <Analysis>
        <Gnuplotter time-step="500" decorate="false">
            <Terminal size="400 400 0" name="png"/>
            <Plot>
                <Field symbol-ref="a" min="0"/>
            </Plot>
        </Gnuplotter>
        <Logger time-step="250">
            <Input>
                <Symbol symbol-ref="a"/>
            </Input>
            <Output>
                <TextOutput file-format="csv"/>
            </Output>
            <Plots>
                <Plot time-step="0">
                    <Style style="lines" line-width="5.0"/>
                    <Terminal terminal="png"/>
                    <X-axis>
                        <Symbol symbol-ref="space.x"/>
                    </X-axis>
                    <Y-axis minimum="0" maximum="3">
                        <Symbol symbol-ref="a"/>
                    </Y-axis>
                    <Color-bar>
                        <Symbol symbol-ref="i"/>
                    </Color-bar>
                    <Range>
                        <Data/>
                        <Time mode="current"/>
                    </Range>
                </Plot>
            </Plots>
            <Restriction>
                <Slice value="size.y/2" axis="y"/>
            </Restriction>
        </Logger>
        <Logger time-step="100">
            <Input>
                <Symbol symbol-ref="a"/>
            </Input>
            <Output>
                <TextOutput file-format="matrix"/>
            </Output>
            <Plots>
                <SurfacePlot time-step="500">
                    <Color-bar>
                        <Symbol symbol-ref="a"/>
                    </Color-bar>
                    <Terminal terminal="png"/>
                </SurfacePlot>
            </Plots>
            <Restriction>
                <Slice value="size.y/2" axis="y"/>
            </Restriction>
        </Logger>
    </Analysis>
</MorpheusModel>


In Morpheus GUI: Examples → PDE → ActivatorInhibitor_2D.xml

Things to try

  • Alter the level of saturation of the activator by changing kappa (for example, using a ParamSweep). By changing this parameter, the model can generate stripes, nets and spots.

Domains: Reaction-diffusion in irregular domains

Spot pattern in Gierer-Meinhardt model with irregular domain

Introduction

A 2D activator-inhibitor model (Meinhardt and Gierer, 1972), solved in a irregular domain that is load from file.

Model description

This model uses an irregular Domain with constant boundary conditions. The domain is loaded from a TIFF image. See Space/Lattice/Domain.

Model

Note: This model requires the external file domain.tif (download domain.tif).

ActivatorInhibitor_Domain.xml

<MorpheusModel version="3">
    <Description>
        <Title>Example-ActivatorInhibitor-2D-Domain</Title>
        <Details>Meinhardt-Gierer (activator-inhibitor) model solved in a nonregular domain with constant boundaries.</Details>
    </Description>
    <Global>
        <Field symbol="a" value="rand_norm(0.5,0.1)" name="activator">
            <Diffusion rate="0.02"/>
            <BoundaryValue boundary="domain" value="0.01"/>
        </Field>
        <Field symbol="i" value="0.1" name="inhibitor">
            <Diffusion rate="0.25"/>
            <BoundaryValue boundary="domain" value="0"/>
        </Field>
        <System solver="runge-kutta" time-step="5" name="Meinhardt">
            <Constant symbol="rho" value="0.001"/>
            <Constant symbol="rho_a" value="0.001"/>
            <Constant symbol="mu_i" value="0.02"/>
            <Constant symbol="mu_a" value="0.04"/>
            <Constant symbol="kappa" value="0.10"/>
            <DiffEqn symbol-ref="a">
                <Expression>(rho/i)*((a^2)/(1 + kappa*a^2)) - mu_a * a + rho_a</Expression>
            </DiffEqn>
            <DiffEqn symbol-ref="i">
                <Expression>rho*((a^2)/(1+kappa*a^2)) - mu_i *i</Expression>
            </DiffEqn>
        </System>
    </Global>
    <Space>
        <Lattice class="square">
            <Size symbol="size" value="100 100 0"/>
            <NodeLength value="1"/>
            <Domain boundary-type="constant">
                <Image path=":/examples/PDE/domain.tif"/>
            </Domain>
            <Neighborhood>
                <Order>1</Order>
            </Neighborhood>
        </Lattice>
        <SpaceSymbol symbol="space"/>
    </Space>
    <Time>
        <StartTime value="0"/>
        <StopTime value="10000"/>
        <SaveInterval value="0"/>
        <RandomSeed value="2"/>
        <TimeSymbol symbol="time"/>
    </Time>
    <Analysis>
        <Gnuplotter time-step="200" decorate="false">
            <Terminal size="400 400 0" name="png"/>
            <Plot>
                <Field symbol-ref="a" min="0"/>
            </Plot>
        </Gnuplotter>
    </Analysis>
</MorpheusModel>



In Morpheus GUI: Examples → PDE → ActivatorInhibitor_Domain.xml


Spatial parameter sweep: Turing patterns

Spots and stripes appear under various conditions in linear Turing model

Introduction

This model shows the pattern formation abilities of Turing's linear activator-inhibitor model (Miyazawa et al., 2010). It shows how to vary parameters as a function of space.

Model description

Instead of fixed parameters defined as Constants, this model uses Functions for two parameters. The parameters C (activator production) and A (rate of auto-activation) as defined as Function of space and varied over the X- and Y-axes respectively. This requires the definition of a SpaceSymbol that can be used in expressions.

The results show the appearance of white spots (left), black spots (right) and labyrinthine patterns (middle).

Model

TuringPatterns.xml

<MorpheusModel version="3">
    <Description>
        <Title>Example-TuringPatterns</Title>
        <Details>
Miyazawa, Okamoto and Kondo, Blending of animal colour patterns by hybridization, Nature Communications, 2010</Details>
    </Description>
    <Global>
        <Field symbol="u" value="4.1+rand_uni(0,1)">
            <Diffusion rate="1"/>
        </Field>
        <Field symbol="v" value="4.84+rand_uni(0,1)">
            <Diffusion rate="20"/>
        </Field>
        <System solver="runge-kutta" time-step="0.25" name="Miyazawa">
            <Function symbol="A">
                <Expression>0.07 + ((0.07 * l.y)/ s.y)</Expression>
            </Function>
            <Constant symbol="B" value="0.08"/>
            <Function symbol="C">
                <Expression>-0.1 + ((0.5 * l.x)/ s.x)</Expression>
            </Function>
            <Constant symbol="D" value="0.03"/>
            <Constant symbol="E" value="0.10"/>
            <Constant symbol="F" value="0.12"/>
            <Constant symbol="G" value="0.06"/>
            <Constant symbol="R" value="20.0"/>
            <Constant symbol="synU_max" value="0.23"/>
            <Constant symbol="synV_max" value="0.50"/>
            <Function symbol="s_u">
                <Expression>max( 0, min( synU_max, A*u-B*v+C))</Expression>
            </Function>
            <Function symbol="s_v">
                <Expression>max( 0, min( synV_max, E*u - F))</Expression>
            </Function>
            <DiffEqn symbol-ref="u">
                <Expression>R*(s_u - D*u)</Expression>
            </DiffEqn>
            <DiffEqn symbol-ref="v">
                <Expression>R*(s_v - G*v)</Expression>
            </DiffEqn>
        </System>
    </Global>
    <Space>
        <Lattice class="square">
            <Size symbol="s" value="512 512 0"/>
            <NodeLength value="1"/>
            <BoundaryConditions>
                <Condition boundary="x" type="noflux"/>
                <Condition boundary="y" type="noflux"/>
            </BoundaryConditions>
            <Neighborhood>
                <Order>1</Order>
            </Neighborhood>
        </Lattice>
        <SpaceSymbol symbol="l"/>
    </Space>
    <Time>
        <StartTime value="0"/>
        <StopTime value="30"/>
        <SaveInterval value="0"/>
        <RandomSeed value="1"/>
        <TimeSymbol symbol="time"/>
    </Time>
    <Analysis>
        <Gnuplotter time-step="2" decorate="false">
            <Terminal persist="true" name="png"/>
            <Plot>
                <Field symbol-ref="u">
                    <ColorMap>
                        <Color value="1" color="black"/>
                        <Color value="0.0" color="white"/>
                    </ColorMap>
                </Field>
            </Plot>
        </Gnuplotter>
        <Logger time-step="0.0">
            <Input>
                <Symbol symbol-ref="u"/>
            </Input>
            <Output>
                <TextOutput file-format="csv"/>
            </Output>
            <Restriction>
                <Slice value="s.x/2" axis="x"/>
            </Restriction>
            <Plots>
                <Plot time-step="-1" name="slice at x is half">
                    <Style style="lines" line-width="3.0"/>
                    <Terminal terminal="png"/>
                    <X-axis>
                        <Symbol symbol-ref="l.y"/>
                    </X-axis>
                    <Y-axis>
                        <Symbol symbol-ref="u"/>
                    </Y-axis>
                    <Color-bar>
                        <Symbol symbol-ref="v"/>
                    </Color-bar>
                </Plot>
            </Plots>
        </Logger>
        <Logger time-step="2">
            <Input>
                <Symbol symbol-ref="u"/>
            </Input>
            <Output>
                <TextOutput file-format="csv"/>
            </Output>
            <Plots>
                <SurfacePlot time-step="2">
                    <Color-bar>
                        <Symbol symbol-ref="u"/>
                    </Color-bar>
                    <Terminal terminal="png"/>
                </SurfacePlot>
            </Plots>
        </Logger>
    </Analysis>
</MorpheusModel>


In Morpheus GUI: Examples → PDE → TuringPatterns.xml


3D reaction-diffusion: Excitable Media

Scroll wave appear in Barkley model of excitable media in 3D

Introduction

This example uses the Barkley model of excitable media, similar to the Fitzhugh-Nagumo model to show how to model and visualize reaction-diffusion models in 3D.

Model description

This model defines a 3D cubic Lattice with noflux BoundaryConditions. Two Layers are defined for the two species: 'u' is the signal, and 'v' the refractoriness. As in the examples above, the DiffEqn as specified in the System in PDE. Nothing strange here.

To visualize the resulting scrolling waves in 3D, the TiffPlotter is used. This Analysis plugin writes TIFF image stacks that cam be opened by image analysis software such as Fiji/ImageJ. To import Morpheus TIFF images into Fiji, macros scripts are available that help you to create 3D (xyz), 4D (xyzt) or even 5D (xyzct) images and movies of your simulations.

Although unable to plot 3D, the GnuPlotter can still be helpful to plot a 2D slice. See Analysis/Gnuplotter/PDE/slice.

Model

ExcitableMedium_3D.xml

<MorpheusModel version="3">
    <Description>
        <Title>Example-ExcitableMedium-3D</Title>
        <Details>Simulates the Barkley model of an excitable medium, see: http://www.scholarpedia.org/article/Barkley_model
Derived from FitzHugh-Nagumo model and Hogdkin-Huxley model.

TIFF images can be viewed with external tools such as Fiji / ImageJ (http://fiji.sc/Fiji) or BioView3D (http://www.dimin.net/software/bioview3d/). The latter also reads the OME header for 3D,4D and 5D images.
 </Details>
    </Description>
    <Global>
        <Field symbol="u" value="if( l.x>=s.x/2-5 and l.x&lt;=s.x/2+5 and l.z>=s.z/2-5 and l.z&lt;=s.z/2+5 and l.y&lt;=s.y/4 , 1, 0 )">
            <Diffusion rate="0.5"/>
        </Field>
        <Field symbol="v" value="if(l.x&lt;=s.x/2 and l.z&lt;=(3*s.z)/4, 1, 0)">
            <Diffusion rate="0.5"/>
        </Field>
        <System solver="runge-kutta" time-step="0.05">
            <DiffEqn symbol-ref="u">
                <Expression>(1/e)*u*(1-u)*(u-((v+b)/a))</Expression>
            </DiffEqn>
            <DiffEqn symbol-ref="v">
                <Expression>u-v</Expression>
            </DiffEqn>
            <Constant symbol="e" value="0.02"/>
            <Constant symbol="a" value="0.8"/>
            <Constant symbol="b" value="0.01"/>
        </System>
    </Global>
    <Space>
        <Lattice class="cubic">
            <Size symbol="s" value="50 50 50"/>
            <BoundaryConditions>
                <Condition boundary="x" type="noflux"/>
                <Condition boundary="y" type="noflux"/>
                <Condition boundary="z" type="noflux"/>
                <Condition boundary="-x" type="noflux"/>
                <Condition boundary="-y" type="noflux"/>
                <Condition boundary="-z" type="noflux"/>
            </BoundaryConditions>
            <NodeLength value="1.0"/>
            <Neighborhood>
                <Order>1</Order>
            </Neighborhood>
        </Lattice>
        <SpaceSymbol symbol="l" name="position in space"/>
    </Space>
    <Time>
        <StartTime value="0"/>
        <StopTime value="25"/>
        <SaveInterval value="0"/>
        <TimeSymbol symbol="time"/>
    </Time>
    <Analysis>
        <Gnuplotter time-step="5">
            <Terminal size="600 600 0" name="png"/>
            <Plot>
                <Field symbol-ref="u" slice="25"/>
            </Plot>
        </Gnuplotter>
        <TiffPlotter timelapse="true" format="32bit" OME-header="true" compression="false" time-step="0.5">
            <Channel symbol-ref="u"/>
            <Channel symbol-ref="v"/>
        </TiffPlotter>
    </Analysis>
</MorpheusModel>


In Morpheus GUI: Examples → PDE → ExcitableMedium_3D.xml

Things to try

  • Import resulting sequence of TIFF images in ImageJ or Fiji, and create 4D movie using ImageJ's 3D plugin:
    1. Open u_v.tif in ImageJ: File → Open.
    2. Create hyperstack: Image → Hyperstack → Convert to Hyperstack. Channels ( c ): 2, Slices (z): 50, Frames (t): 51, Display Mode: Composite.
    3. Display in 4D: Plugins → 3D Viewer. Use default parameters. Press OK.

examples/reaction-diffusion.txt · Last modified: 14:17 02.10.2013 by Walter