Solving the Gierer-Meinhardt model using Julia

The Gierer-Meinhardt model is defined as follows:

$$ \frac{\partial u}{\partial t} = D_u \Delta u + \rho\frac{u^2}{v} - \mu_u u + \rho_u $$ $$ \frac{\partial v}{\partial t} = D_v \Delta v + \rho u^2 - \mu_v v + \rho_v, $$

with

$u$ being a short-range autocatalytic substance, or in other words, an activator,
$v$ being its long-range antagonist, or in other words, an inhibitor, and,
$\Delta = \sum\limits_{i = 1}^{n} \frac{\partial^2}{\partial x_i^2}$ being the n-dimensional Laplace operator.

In the Gierer-Meinhardt equations, the autocatalytic substance activates both itself and the inhibitor substance (with rate $\rho u^2)$, whereas the inhibitor function inhibits the growth of the autocatalytic substance (with rate $\frac{1}{v})$. Both substances have a natural decay rate of the form $\mu_u u$ and $\mu_v v$ respectively. Finally, both substances have an activator-independent production rate ($\rho_u$ and $\rho_v$).

For the right choice of parameters, pattern formation can be observed in the solution of the Gierer-Meinhardt model.

More info on the Gierer-Meinhardt model can be found in this scholarpedia article.

Solving the system

We’ll solve the following IBVP:

$$ \frac{\partial u}{\partial t} = D_u \Delta u + \rho\frac{u^2}{v} - \mu_u u + \rho_u $$ $$ \frac{\partial v}{\partial t} = D_v \Delta v + \rho u^2 - \mu_v v + \rho_v $$ $$ u(x,y,0) = \exp{\left(-(x-a)^2-(y-a)^2\right)}, \quad \forall (x,y) \in (0,L)^2$$ $$ v(x,y,0) = {\rm rand}(), \quad \forall (x,y) \in (0,L)^2$$ $$ \frac{\partial u}{\partial n}(x,0,t) = \frac{\partial v}{\partial n}(x,0,t) = 0, \quad \forall x \in [0,L]$$ $$ \frac{\partial u}{\partial n}(x,L,t) = \frac{\partial v}{\partial n}(x,L,t) = 0, \quad \forall x \in [0,L]$$ $$ \frac{\partial u}{\partial n}(0,y,t) = \frac{\partial v}{\partial n}(0,y,t) = 0, \quad \forall y \in [0,L]$$ $$ \frac{\partial u}{\partial n}(L,y,t) = \frac{\partial v}{\partial n}(L,y,t) = 0, \quad \forall y \in [0,L] .$$

We start by initializing the parameters

Du = 1;
Dv = 100;

ρ_u = 0.5;
ρ_v = 0;
ρ = 1;
μ_u = 1;
μ_v = 6.1;

a = 5;

L = 100;

Then, we discretize the 2D space grid

sizez = L;  # size of the 2D grid
dx = 100.0 / sizez;  # space step

and then we discretize time

T = 150.0;  # total time
dt = 0.0003;  # time step
n = floor(Int, (T / dt));  # number of iterations

nvis = 500; # saves every nvis time steps

Next we use a random seed in order for the results to be the same across multiple runs and set the initial conditions

using Random
Random.seed!(1)
U = [exp(-(x-a)^2-(y-a)^2) for x in 1:1:sizez, y in 1:1:sizez];
V = rand(sizez, sizez);

Consequently, we initialize two matrices for the inside of the grid, and five matrices that will help us compute the Laplacian

Uc = zeros(Float64, sizez - 2, sizez - 2);
Vc = copy(Uc);
Ztop = copy(Uc);
Zleft = copy(Uc);
Zbottom = copy(Uc);
Zright = copy(Uc);
Zcenter = copy(Uc);

Finally, we initialize a 3D matrix, in order to save $u$ every $nvis$ time steps

UR = zeros(sizez, sizez, floor(Int, n/nvis));

We define a function in order to discretize the Laplacian, using the broadcasting abilities of Julia and the @views macro.

@views function DLaplacian(Z) # Centered differences discretization of the Laplacian
    Ztop .= Z[1:end-2, 2:end-1]
    Zleft .= Z[2:end-1, 1:end-2]
    Zbottom .= Z[3:end, 2:end-1]
    Zright .= Z[2:end-1, 3:end]
    Zcenter .= Z[2:end-1, 2:end-1]
    return (Ztop .+ Zleft .+ Zbottom .+ Zright .- 4 .* Zcenter) ./ (dx^2)
end

DLaplacian (generic function with 1 method)

We then proceed to iterate the solution beginning from the second timestep and ending in timestep $n$.

@views begin
    j = 1
    for i in 2:n
        Uc .= U[2:end-1, 2:end-1]
        Vc .= V[2:end-1, 2:end-1]
    
        U[2:end-1, 2:end-1] .= Uc .+ dt .* (Du .* DLaplacian(U) .+ ρ_u .- μ_u .* Uc .+ ρ.* (Uc.^2)./Vc )
        V[2:end-1, 2:end-1] .= Vc .+ dt .* (Dv .* DLaplacian(V) .+ ρ_v .- μ_v .* Vc .+ ρ.* Uc.^2 )
    
        for Z in (U, V)  # Neumann boundary conditions
            Z[1, :] .= Z[2, :]
            Z[end, :] .= Z[end-1, :]
            Z[:, 1] .= Z[:, 2]
            Z[:, end] .= Z[:, end-1]
        end

        if i%nvis == 0
            UR[:,:,j] .= U
            j += 1
        end
    end
end

We notice interesting spot-like patterns in the solution for $t = 1000 \cdot dt \cdot nvis$

using CairoMakie

joint_limits = (0, 17)

time_step = 1000
time = time_step * dt * nvis

fig = Figure(backgroundcolor = "#9B89B3",
             resolution = (600, 600))
ax = Axis(fig[1,1],
          title = "Solution of the Gierer-Meinhardt model \n for t=$(time)",
          xlabel = "x",
          ylabel = "y")
hmr = CairoMakie.heatmap!(ax, UR[:,:,time_step])
Colorbar(fig[:,end+1], colorrange = joint_limits, label = L"U")
fig

Of course, if we want to speed things up, we’ll want to wrap the above code in a function

using Random

@views function GiererMeinhardt() 
    Du = 1;
    Dv = 100;
    
    ρ_u = 0.5;
    ρ_v = 0;
    ρ = 1;
    μ_u = 1;
    μ_v = 6.1;
    
    a = 5;
    
    L = 100;

    nvis = 1000
    
    sizez = 100  # size of the 2D grid
    dx = 100.0 / sizez  # space step
    
    T = 150.0  # total time
    dt = 0.0003  # time step
    n = floor(Int, (T / dt))  # number of iterations
    
    Uc = zeros(Float64, sizez - 2, sizez - 2)
    Vc = copy(Uc)
    Ztop = copy(Uc)
    Zleft = copy(Uc)
    Zbottom = copy(Uc)
    Zright = copy(Uc)
    Zcenter = copy(Uc)

    Random.seed!(1)
    U = [exp(-(x-a)^2-(y-a)^2) for x in 1:1:sizez, y in 1:1:sizez]
    V = fill(1.0, sizez, sizez)

    UR = zeros(sizez, sizez, floor(Int, n/nvis))
    
    function DLaplacian(Z) # Centered differences discretization of the Laplacian
        Ztop .= Z[1:end-2, 2:end-1]
        Zleft .= Z[2:end-1, 1:end-2]
        Zbottom .= Z[3:end, 2:end-1]
        Zright .= Z[2:end-1, 3:end]
        Zcenter .= Z[2:end-1, 2:end-1]
        return (Ztop .+ Zleft .+ Zbottom .+ Zright .- 4 .* Zcenter) ./ (dx^2)
    end
    
    j = 1
    for i in 2:n
        Uc .= U[2:end-1, 2:end-1]
        Vc .= V[2:end-1, 2:end-1]
    
        U[2:end-1, 2:end-1] .= Uc .+ dt .* (Du .* DLaplacian(U) .+ ρ_u .- μ_u .* Uc .+ ρ.* (Uc.^2)./Vc )
        V[2:end-1, 2:end-1] .= Vc .+ dt .* (Dv .* DLaplacian(V) .+ ρ_v .- μ_v .* Vc .+ ρ.* Uc.^2 )
    
        for Z in (U, V)
            Z[1, :] .= Z[2, :]
            Z[end, :] .= Z[end-1, :]
            Z[:, 1] .= Z[:, 2]
            Z[:, end] .= Z[:, end-1]
        end

        if i%nvis == 0
            UR[:,:,j] .= U
            j += 1
        end
    end
    return UR
end

resG = GiererMeinhardt();

Finally, we make a nice little video illustrating how the pattern formed

using CairoMakie

joint_limits = (0, 17)

fig = Figure(backgroundcolor = "#9B89B3",
             resolution = (600, 600))
ax = Axis(fig[1,1],
          title = "Solution of the Gierer-Meinhardt model",
          xlabel = "x",
          ylabel = "y")
hmr = CairoMakie.heatmap!(ax, resG[:,:,1])
Colorbar(fig[:,end+1], colorrange = joint_limits, label = L"U")
fig

nframes = 10
framerate = 30
iterator = 1:2:500

output_gif = record(fig, "GiererMeinhardt.mp4", iterator;
        framerate = framerate) do t
    CairoMakie.heatmap!(ax, resG[:,:,t], colorrange = joint_limits)
end

Of course, different initial conditions and parameter values, will change the behavior of the model. So, copy the code and do some experimenting!