# Brusselator and Catalyst

Good evening to all,

I recently solved the ODEs system for the brusselator:

``````A -> X (k1)
2X + Y -> 3X (k2)
B + X -> Y + D (k3)
X -> E (k4)
``````

where A, B, etc are chemical species and k represents kinetics coefficients.
Because I wanted to get the limit cycle ([Y] as a function of [ X]), I chose:

``````k1 = k2 = k3 = k4 = 1
[X]0 = [Y]0 = 1
[A] = 1 and [B] = 3
``````

Here [ ] denotes the concentration in various species. The objective was to solve d[ X]/dt and d[Y]/dt. It is not difficult to fix the value of [A] and [B] since there is a direct access for the ODEs to be modified:

``````d[X]/dt = [A] + [Y][X]^2 - [B][X] -[X]
d[Y]/dt = -Y][X]^2 + [B][X]
``````

I then became interested in Catalyst (so I solved the thermal decomposition of ethane as well as a problem of 2 consecutive organic reactions).
The next problem is to solve the brusselator with Catalyst. I defined:

``````brusselator = @reaction_network begin
k1, A → X
k2, 2X + Y → 3X
k3, B + X → Y + D
k4, X → E
end k1 k2 k3 k4

``````

I disabled rescaling of reaction rates. The latexification shows that the ODEs so obtained are correct.
There is no problem to write p and tspan.

However how to formulate u0 since [A] and [B] are not initial concentrations? How to introduce “external” modifications to the ODEs system since this latter is no longer apparent?
At first sight I thought that some kind of forcing could help me. However I don’t see how to proceed on the basis of examples given here.

Thierry

If `A` and `B` are constant, those reactions are really zero and first order respectively. So you can just move `A` and `B` into the rates like

``````brusselator = @reaction_network begin
k1*A, 0 → X
k2, 2X + Y → 3X
k3*B, X → Y + D
k4, X → E
end k1 k2 k3 k4 A B
``````

Alternatively, you need to modify your reactions to ensure they don’t change `A` and `B`, then the ODEs for `A` and `B` will have zero derivatives so you can just set their values in the initial condition vector:

``````brusselator = @reaction_network begin
k1, A → A + X
k2, 2X + Y → 3X
k3, B + X → B + Y + D
k4, X → E
end k1 k2 k3 k4
``````

edit: Notice, in the first example `A` and `B` are now parameters.

2 Likes

Hello Sam,

``````brusselator = @reaction_network begin
k1*A, 0 → X
k2, 2X + Y → 3X
k3*B, X → Y + D
k4, X → E
end k1 k2 k3 k4 A B
``````

with p = (1.0,1.0,1.0,1.0,1.0,3.0)
u0 = [1.0,1.0,0.0,0.0]
tspan=(.0,50.0)

perfectly works.
Many thanks!  [I’m sorry I didn’t think to rework the system of equations.]  2 Likes