Solving the Gierer-Meinhardt model using Julia
Learning how to discretize and solve a reaction-diffusion system in Julia
![](/post/gierermeinhardt/featured_hu20ece05c6429352de109630aa7607bd5_47629_720x2500_fit_q75_h2_lanczos_3.webp)
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, download the notebook and do some experimenting!