Brusselator and Catalyst

Specific Domains Chemistry
1

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.

Could someone please guide me?
Thank you in advance,
Thierry

2

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
3

Hello Sam,

Your solution:

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! :wink: :+1:
[I’m sorry I didn’t think to rework the system of equations.]

brusselator_1 brusselator_2

2 Likes
4

Glad that helped!