Brusselator

From bio-physics-wiki

Jump to: navigation, search
Boris Pavlovich Belousov. He studied chemistry at the ETH Zurich, but got no degree for financial resons. He continued doing science back in Russia where he discovered his famous reaction.
The Brusselator is a system of nonlinear ODEs, that describe the dynamics of the Belousov-Zhabotinsky reaction. While working on the citric acid cycle, Belousov studied a mixture of citric acid and bromate ions in a solution of sulfuric acid in presence of a cerium catalyst. Along with the Briggs-Rauscher reaction (~1921), this was one of the first reactions, where biochemical oscillations have been observed (~1958). At this time the view point, that entropy has to increase continuously by the second law of thermodynamics until entropy reaches its maximum, was so dominant, that the observations which Belousov made, were thought to be an artefact. Thus Belousov couldn't manage to get his work published, till 1959 when it came out in the obscure proceedings of a Russian medical meeting (from Strogatz). In 1961 Zhabotinsky reproduced Belousovs work showed further oscillating reactions. Only in 1968 the results were shown to the western scientific community at a conference in Prague, since the political climate was not the best. The view point that entropie has to increase monotonically until a maximum is reached was shown to be restricted to closed Systems. Open systems instead can give rise to coherent behaviour such as macroscopic flows (Bénard instability), space-time patterns and oscillations. Therefore far-from-equilibrium conditions can be a source of order and lead to self-organisation. Since the maintainance of far-from-equilibrium conditions require a continuous flow of engergy such systems are also calles dissipative-systems. The first theoretical explanation of oscillating chemical reactions was accomplished by Ilya Prigogine and René Lefever in 1968 in the trimolecular model also called the Brusselator.

Nicolis & Prigogine

Firstly, it is important to realize that the primary objective in studying a model is to discover the types of qualitative behavior compatible with some fundamental laws, such as the laws of thermodynamics and chemical kinetics. In this respect, regardless of the reasonableness or not of a tri- molecular reaction step, the Brusselator is a perfectly acceptable model for the study of cooperative processes in chemical kinetics. It plays somewhat the same role as such models as the harmonic oscillator, or the Heisenberg model in ferromagnetism, which have both been widely used to illustrate basic features of systems described by the laws of classical or quantum mechanics. (in Self-Organization in Nonequilibrium Systems: From Dissipative Structures to Order through Fluctuations)

The reaction mechanism for the Brusselator \begin{equation} A \overset{k_1}{\underset{k_{-1}}{\rightleftharpoons}} X\\ B+X \overset{k_2}{\underset{k_{-2}}{\rightleftharpoons}} Y+D\\ 2X+Y \overset{k_3}{\underset{k_{-3}}{\rightleftharpoons}} 3X\\ X \overset{k_4}{\underset{k_{-4}}{\rightleftharpoons}} E\\ \end{equation} If we assume that no backward reactions occure and that the final products $D,E$ are removed immediatley after producion, this means $k_{-2}=0$ and $k_{-4}=0$ and symlifying the system further by assuming $k_{-1}=k_{-3}=0$ the dynamical system reads \begin{align} \frac{dX}{dt}&=1-(b+1)X+aX^2Y\\ \frac{dY}{dt}&=bX-aX^2Y \tag{1} \end{align} We obtain the fixed point by setting the system of simultaneous nonlinear differential equations (1) to zero. \begin{align} 1-(b+1)X+aX^2Y&=0\\ bX-aX^2Y&=0 \tag{2} \end{align} Adding the equations of (2) yields $x=1$. Plugging this result into one of the equations of (2) gives $y=b/a$. So there is only one singular point \begin{align} (x^*,y^*)=\left(1,\frac{b}{a}\right) \end{align} What about the stability of this fixed point? To answer this we compute the Jacobian \begin{align} Df(x,y)=\begin{pmatrix}1-b + 2aXY & aX^2\\ b-2aXY &aX^2 \end{pmatrix} \end{align} and plug in the values for the fixed point. \begin{align} Df(x,y)=\begin{pmatrix}1-b & a\\ -b &a \end{pmatrix} \end{align} To determine the solution of the charakteristic polynomial we compute the trace $\tau$ and the determinant $\Delta$ of the Jacobian \begin{align} \tau &=trace(Df(1,b/a))=1-b+a\\ \Delta&=det(Df(1,b/a))=a \end{align} Now notice, for $\tau >0$ the fixed point is unstable, while for $\tau <0$ it is stable, changing the stability at $\tau =0$. So we find the critical parameter values by the equation \begin{align} \tau &=1-b+a=0\\ \end{align}

Parameterspace for the Brusselator

If you plug in the point $(a,b)=(0,0)$ of the parameterspace, then $\tau$ becomes positive and is therefore unstable at this point. Hence, the region under the curve $b=1+a$ is unstable. To sketch the phaseportait we first find the nullclines \begin{align} 1-(b+1)X+aX^2Y&=0\\ \Rightarrow Y=\frac{(b+1)X-1}{aX^2} \end{align} \begin{align} bX-aX^2Y&=0 \\ \Rightarrow Y=\frac{b}{aX} \end{align}

Brusselatoranalyze.png

To proof the existence of a closed orbit we use the Poincaré-Bendixon-Theorem and construct a trapping region. First the condition $\tau<0$ must be satisfied, then the fixed point $(1,b/a)$ is unstable and we can find a region around the fixed point where all trajecotries point away from this region.

Trapping1.png

We use the x-axis as boundery for the trapping region and start another line paralell to y-axis, starting from the intersection of the nullcline with the x-axis. Then we can be sure that all vectors on the boundery point to the right. \begin{align} \frac{(b+1)X-1}{aX^2}=0\\ \Rightarrow X=\frac{1}{1+b} \end{align} The boudery region on "the top" can be found at the intersection of $X=\frac{1}{1+b}$ with the nullcline from $\dot{y}=0$ and is parallel to the x-axis. \begin{align} Y&=\frac{b}{aX} \\ \Rightarrow Y&=\frac{b(b+1)}{a} \end{align} Now we chose the right boundery for $x=k>1$ which reaches up till it intersects with the nullcline from $\dot{x}=0$. We have to choose $k$ large enough so that the vectors on the diagonal, joining the top and right boundery, pointing in the trapping region. This is the case, when the product of the normal vector $\mathbf{n}=(n_1,n_2)$ with the vectorfield (1) is negativ. \begin{align} \begin{pmatrix} n_1 \\ n_2 \end{pmatrix} \cdot \begin{pmatrix} \dot x \\ \dot y \end{pmatrix} < 0 \end{align} We simply choose $n_1=n_2=1$ and get \begin{align} 1-(b+1)X+aX^2Y+bX-aX^2Y<0\\ 1<X \end{align} So as we choose $k$ large enough, the condition $X>1$ is fullfilled anywhere on the dashed line.

Trapping2.png

Since all vectors on the boundery point into the trapping region there must exist a stable limit cycle by the Poincaré-Bendixon-Theorem. Finally we want to know the time period $T=2\pi/\omega$ a trajectory needs for one oscillation near the critical parameter $b_c=a+1$ where $\tau=0$. To this end we calculate the imaginary part of $\lambda$. So we get for the angular frequency $\omega$ \begin{align} i \omega=|Im(\lambda_{c})|= \left| Im \left(\frac{\tau \pm \sqrt{\tau^2-4 \Delta}}{2} \right)\right|=\frac{\sqrt{-4 \Delta}}{2}=i\sqrt{a} \end{align} this means the trajectory runs through one period in \begin{align} T=\frac{2 \pi}{\sqrt{a}} \end{align}