Bachelor’s Thesis: The Brusselator

“Development of a Python program for investigation of reaction-diffusion coupling in oscillating reactions”

Github

The Brusselator is a theoretical model for a type of autocatalytic reaction. This project aims aims to solve the brusselator equations and apply them to a two-dimensional cellular automaton.
https://github.com/lukaswittmann/brusselator
0 forks.
1 stars.
0 open issues.

Recent commits:

Add files via uploadtypes and utils by certik (https://github.com/certik/fortran-utils.git), GitHub
Create main.f90, GitHub
Delete constants.f90, GitHub
Delete main.f90, GitHub
Delete plot.py, GitHub

Introduction

The Brusselator is a model developed by I. Prigogine and R. Lefever at the Universite Libre de Bruxelles in Belgium. [16] developed for the description of nonlinear systems. The model includes some assumptions [12, 11], such as irreversible reactions ( $k_{-i}=0$ for $i=1,2,3,4)$ and constant reactant and product species concentrations A, B, C and D. Realistically, this happens by constant addition of A and B and removal of C and D.

$\ce{A} \ce{->[k_1]} \ce{X}$

$\ce{B + X} \ce{->[k_2]} \ce{Y + D}$

$\ce{2 X + Y} \ce{->[k_3]} \ce{3 X}$

$\ce{X} \ce{->[k_4]} \ce{E}$

$\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_$

$\ce{A + B} \ce{->} \ce{D + E}$

With the assumptions of chapter 2.1, the dierential kinetic time laws of the species X and Y can be derived. The concentration of a species is noted with square brackets. The concentration notation is assumed to be dimensionless in the following, since there is a mathematical consideration of the issue. With $[D]=[E]=0$ and $k_{-i}=0$ for $i=1,2,3,4$ the differential time laws follow.

$f\left([X],[Y]\right) \;\equiv\; \dfrac{\text{d} \left[X\right]}{\text{d} t} = k_1 \left[A\right] - k_2 \left[B\right] \left[X\right] + k_3 \left[X\right]^2 \left[Y\right] - k_4 \left[X\right] \label{eq_DGL_f}$

$g\left([X],[Y]\right) \;\equiv\; \dfrac{\text{d} \left[Y\right] }{\text{d} t} = k_2 \left[B\right] \left[X\right] - k_3 \left[X\right]^2 \left[Y\right] \label{eq_DGL_g}$

Spatially Homogeneous Brusselator

Dynamic Analysis of the Spatially Homogeneous Brusselator

This part is rather long and thus shortened. The full version is shown in detail in my thesis. The fixpoint can be calculated.

$\left[X \right]^* = \dfrac{k_1 \left[A\right]}{k_4} \label{eq_fixpunkt_X}$

$\left[Y \right]^* = \dfrac{k_2 k_4 \left[B\right]}{k_3 k_1 \left[A\right]}$

$F=\left(\dfrac{k_1 \left[A\right]}{k_4} ,\; \dfrac{k_2 k_4 \left[B\right]}{k_3 k_1 \left[A\right]}\right)$

It is not possibile to analyse non-linear differential equations. Because of this, the System needs to be linearised by using an approximation. The simplified Idea behind my procedure was using a disturbance and investigating its behaviour.

$U = [X] - [X]^*$

$V = [Y] - [Y]^*$

$\dfrac{\text{d} U}{\text{d} t}= \dfrac{\text{d} [X]}{\text{d} t}$

$\dfrac{\text{d} V}{\text{d} t} = \dfrac{\text{d} [Y]}{\text{d} t}$

Next both functions $f\left([X],[Y]\right)$ and $g\left([X],[Y]\right)$ are substituted and a Taylor series approximation is used.

$\dfrac{\text{d} U}{\text{d} t} = f\left([X]^* + U, [Y]^* + V\right) = f([X]^*,[Y]^*) + U\dfrac{\delta f([X]^*,[Y]^*)}{\delta [X]}+ V\dfrac{\delta f([X]^*,[Y]^*)}{\delta [Y]} + h\left(U^2,V^2,UV\right)$

Same for $V$ .

$\dfrac{\text{d} U}{\text{d} t} = U\dfrac{\delta f([X]^*,[Y]^*)}{\delta [X]}+ V\dfrac{\delta f([X]^*,[Y]^*)}{\delta [Y]} + h\left(U^2, V^2, UV\right)$

$\dfrac{\text{d} V}{\text{d} t} = U\dfrac{\delta g([X]^*,[Y]^*)}{\delta [X]}+ V\dfrac{\delta g([X]^*,[Y]^*)}{\delta [Y]} + h\left(U^2, V^2, UV\right)$

By neglecting any higher therms ( $h\left(U^2, V^2, UV\right)$ ) it is obvious, that the matrix of partial derivatives is the Jacobian. This is only possible as long as the position and type of the fixed point of the linearized system is also exactly the one of the non-linear system.

$\renewcommand{\arraystretch}{2} \begin{bmatrix} \dfrac{\text{d} U}{\text{d} t} \\ \dfrac{\text{d} V}{\text{d} t} \end{bmatrix} = \begin{bmatrix} \dfrac{\delta f}{\delta \left[X\right]} \dfrac{\delta f}{\delta \left[Y\right]} \\ \dfrac{\delta g}{\delta \left[X\right]} \dfrac{\delta g}{\delta \left[Y\right]} \end{bmatrix} \begin{bmatrix} U \\ V \end{bmatrix}= J_{g,f}\left(\left[X\right],\left[Y\right] \right) \begin{bmatrix} U \\ V \end{bmatrix} = \begin{bmatrix} - k_2 \left[B\right] + 2 k_3 \left[X\right]\left[Y\right] - k_4 & k_3 \left[X\right]^2 \\ k_2 \left[B\right] - 2 k_3 \left[X\right]\left[Y\right] & - k_3 \left[X\right]^2 \end{bmatrix}\begin{bmatrix} U \\ V \end{bmatrix}$

Eigenvalues can easily be calculated and simplified using trace and determinant.

$\lambda_{1,2}=\left[ \begin {array}{c} {\dfrac {-{[A]}^{2}+[B]-1}{2}}+{\dfrac {1}{2}\sqrt {{[A]}^{4}-2\,{[A]}^{2}[B]-2\,{[A]}^{2}+{[B]}^{2}-2\,[B]+1}}\\ \dfrac {-{[A]}^{2}+[B]-1}{2}+{\dfrac {1}{2}\sqrt {{[A]}^{4}-2\,{[A]}^{2}[B]-2\,{[A]}^{2}+{[B]}^{2}-2\,[B]+1}}\end {array} \right]= \dfrac{\tau \pm \sqrt{\tau^2-4 \Delta}}{2}$

Now, the stability and the type of the fixpoint can be calculated. For the spatially homogeneous Brusselator there are four domains depending on [A] and [B]. The domains are shown in Tab. 2.1 and Fig. 2.6. For the spatially homogeneous Brusselator there are four domains depending on [A] and [B]. The domains are shown in Tab. 2.1 and Fig. 2.6. The system behavior is independent of starting conditions $[X]_0$ and $[Y]_0$ as long as $[X]_0\neq[X]^*$ and $[Y]_0\neq[Y]^*$ holds. For $[X]_0=[X]^*$ and $[Y]_0=[Y]^*$ the system would already be in equilibrium.

Figure 2.6: [A]-[B] phase space plot with marked regions providing stable fixed points and unstable fixed points as well as with the marked region resulting in complex eigenvalues.

Limit Cycle

With the use of Poincaré–Bendixson theorem, it can be shown, that a limit cycle is possibile. This depends on the domain present. The proof for the Brusselator model is carried out in [25, p. 42] or more generally in [26, p. 196].

Figure 2.7: Poincare-Bendixon theorem schematically applied to the phase space plot Chap. 5.1 with constraints [A] = 1 and [B] = 3.

Reaction-Diffusion System of Brusselator

Both differential equations can now be advanced by using the second Fick’s law.

$\dfrac{\text{d} \left[X\right]}{\text{d} t} = k_1 \left[A\right] - k_2 \left[B\right] \left[X\right] + k_3 \left[X\right]^2 \left[Y\right] - k_4 \left[X\right] + D_x \nabla^2 \left[X\right]$

$\dfrac{\text{d} \left[Y\right] }{\text{d} t} = k_2 \left[B\right] \left[X\right] - k_3 \left[X\right]^2 \left[Y\right] + D_y \nabla^2 \left[Y\right]$

There is no analytical solution for the RDS of Brusselator.

Dynamic Analysis of the Reaction-Diffusion System

It is found, that the Fixpoint equals the one of the spartially homogenous system.

$F=\left(\dfrac{k_1 \left[A\right]}{k_4} ,\; \dfrac{k_2 k_4 \left[B\right]}{k_3 k_1 \left[A\right]}\right)$

The linearisation works the same way as with the spartially homogenous system, just a little more complicated.

$\dfrac{\text{d} \left[X\right]}{\text{d} t} = k_1 \left[A\right] - k_2 \left[B\right] \left[X\right] + k_3 \left[X\right]^2 \left[Y\right] - k_4 \left[X\right] + D_x \nabla^2 \left[X\right] \quad\equiv\quad f\left([X],[Y]\right) + D_x \nabla^2 \left[X\right]$

$\dfrac{\text{d} \left[Y\right] }{\text{d} t} = k_2 \left[B\right] \left[X\right] - k_3 \left[X\right]^2 \left[Y\right] + D_y \nabla^2 \left[Y\right] \quad\equiv\quad g\left([X],[Y]\right) + D_y \nabla^2 \left[Y\right]$

A disturbance as follows is beeing used.

$[X](r,t) \: \Rightarrow \: [X]^*+\Delta [X](r,t)\: =\: [X]^* + \sin(\omega r + \phi)\Delta [X](t)$

With the use of $S=\sin(\omega r + \phi)$ , $\dfrac{\delta [X]^*}{\delta t}=0$ and $\dfrac{\delta^2 S}{\delta r^2}=-\omega^2S$

$S\dfrac{\delta\Delta[X]}{\delta t}= f\left([X]^*+S\Delta[X],[Y]^*+S\Delta[Y]\right) - D_x \omega^2 S \Delta[X]$

and

$S\dfrac{\delta\Delta[Y]}{\delta t}= g\left([X]^*+S\Delta[X],[Y]^*+S\Delta[Y]\right) - D_y \omega^2 S \Delta[Y]$

will be the result. With the use of $\renewcommand{\arraystretch}{1.2}[C]=\begin{bmatrix} [X] \\ [Y] \end{bmatrix}$ , $$\renewcommand{\arraystretch}{1.2}k\left([C]^*+\sin(\omega r + \phi)\Delta[C]\right)=\begin{bmatrix}f\left([X]^*+S\Delta[X],[Y]^*+S\Delta[Y]\right) \\ g\left([X]^*+S\Delta[X],[Y]^*+S\Delta[Y]\right) \end{bmatrix}$ and $D=\begin{bmatrix} D_x & 0 \\ 0 & D_y \end{bmatrix}$

$\sin(\omega r + \phi)\dfrac{\delta \Delta [C]}{\delta t}= k\left([C]^*+\sin(\omega r + \phi)\Delta[C]\right)-D\omega^2\sin(\omega r + \phi) \Delta[C]$

can be formed. Because all reaction terms are spratially independent and without any spartial operators the vector functin can be aproximated as the Jacobian.

$\renewcommand{\arraystretch}{2}k\left([C]^*+\sin(\omega r + \phi)\Delta[C]\right) \approx k\left([C]^*\right)+\begin{bmatrix}\dfrac{\delta f([X],[Y])}{\delta [X]} & \dfrac{\delta f([X],[Y])}{\delta [Y]}\\\dfrac{\delta g([X],[Y])}{\delta [X]} & \dfrac{\delta g([X],[Y])}{\delta [Y]} \end{bmatrix}\sin(\omega r + \phi)\Delta[C]$

This can be shortened.

$\sin(\omega r + \phi)\dfrac{\delta \Delta [C]}{\delta t}= \sin(\omega r + \phi) J_{g,f}\left(\left[X\right],\left[Y\right] \right)\Delta[C] - D\omega^2\sin(\omega r + \phi) \Delta[C]$

$\dfrac{\delta \Delta [C]}{\delta t} = (J_{g,f}\left(\left[X\right],\left[Y\right] \right)-D\omega^2)\Delta[C]$

The Stability of $[C]^*$ can now be aproximated with the use of

$J_{g,f}\left(\left[X\right],\left[Y\right] \right)-D\omega^2 =\renewcommand{\arraystretch}{2}\begin{bmatrix}\dfrac{\delta f([X],[Y])}{\delta [X]} & \dfrac{\delta f([X],[Y])}{\delta [Y]}\\\dfrac{\delta g([X],[Y])}{\delta [X]} & \dfrac{\delta g([X],[Y])}{\delta [Y]}\end{bmatrix} - \begin{bmatrix} D_x\omega^2 & 0 \\ 0 & D_y\omega^2 \end{bmatrix} \text{.}$

The domains of the spartially homogenous Brusselator are now extended. Four different functions can be found, wich sepreate the domains. This time, two of them refer to the $D_x$ – $D_y$ -phase space.

$D_{y}^{grenz,1}\left([A],[B],D_x \right) = \dfrac{\left([B]+1\pm 2\sqrt{[B]}\right)[A]^2D_x}{[B]^2-2[B]+1}$