Schlögl model - Stationary States far from Equilibrium

Question

The Schlögl model can be represented as follows:

$$ A + 2X \underset{k_2}{\stackrel{k_1}{\rightleftharpoons}} 3X \, , \\ X \underset{k_4}{\stackrel{k_3}{\rightleftharpoons}} B \, . $$

The chemical reactions can be written in ODE form:

$$ \dot{C}_{X} = k_{1}C^{0}_{A}(C_{X})^{2} - k_{2}(C_{X})^{3} - k_{3}C_{X} + k_{4}C_{B}^{0} \, , $$ where we fixed the concentrations of $A$ and $B$ to the initial ones $\left[(C_{A} \equiv C_{A}^{0}) ; (C_{B} \equiv C_{B}^{0})\right]$.

The concentration of stationary states does not evolve, therefore $\dot{C}_{X} = 0$:

$$ k_{1}C^{0}_{A}(C_{X})^{2} - k_{2}(C_{X})^{3} - k_{3}C_{X} + k_{4}C_{B}^{0} = 0 \\ \iff (C_{X})^{3} - \frac{k_{1}}{k_{2}}C^{0}_{A}(C_{X})^{2} + \frac{k_{3}}{k_{2}}C_{X} - \frac{k_{4}}{k_{2}}C_{B}^{0} = 0 \\ \iff k_{2} = \frac{a}{(C_{X})} - \frac{k_{3}}{(C_{X})^{2}} + \frac{b}{(C_{X})^{3}} \, $$ where $a = k_{1}C_{A}^{0}$ and $b = k_{4}C_{B}^{0}$.

In order to find the stationary states, we proceed the following manner:

$$ \frac{dk_{2}}{dC_{X}} = 0 \iff C_{x} = \frac{1}{a}\left(k_{3} \pm \sqrt{k_{3}^{2} - 4ab}\right) $$

If $k_{3}^{2} > 4ab \, \, (1)$, there is 3 possible values and if $k_{3}^{2} < 4ab \, \,(2)$ there is one possible value.

Reichl's book, "A Modern Course in Statistical Mechanics", says that the Schlögl model, away from equilibrium, has multiple stationary states, which seem to represent situation (1) and in equilibrium seem to be (2). However, why does (1) represent the model far from equilibrium and (2) the model in equilibrium?

Answer 1

Ilie et al give a numerical example that has two steady states. Here are the parameters (ignore the propensities, I think they are for the stochastic models):

with $C_A$ set to $\pu{1E5}$ and $C_B$ set to $\pu{2E5}$. Depending on the value of the initial concentration of intermediate, you get one of two possible steady states (figure from the same paper, I did not numerically solve myself):

In our context, the y-axis shows the concentration of the intermediate X. If you want to know whether the system is capable of reaching equilibrium (with fixed concentrations of A and B), use the criterion

$$ \frac{C_B}{C_A} \stackrel{?}{=} K = \frac{k_1 k_3}{k_2 k_4}\tag{3}$$

For the given set of parameters, this is not the case, so it is not at equilibrium (also, you would not expect two possible steady state concentrations of X if it were at equilibrium).

[OP] ... there is 3 possible values...

That seems too much for a quadratic equation. According to Vellela and Qian (2008, doi:10.1098/rsif.2008.0476),

Because the ODE form of Schlögl's model is a cubic, there can be one, two or three steady states for a given set of parameters.

The OP does start out with an expression containing concentrations cubed, so unless some terms cancel out during differentiation, I would expect a different sort of solution for the steady state concentration.

[OP] However, why does (1) represent the model far from equilibrium and (2) the model in equilibrium?

You can check by testing one of the following:

$$ \frac{C_X}{C_A} \stackrel{?}{=} \frac{k_1}{k_2}\tag{4}$$

or

$$ \frac{C_B}{C_X} \stackrel{?}{=} \frac{k_3}{k_4}\tag{5}$$

as long as (3) is also true. I'm not attempting this because I have trouble understanding the derivation and the result given by the OP. Specifically, I am puzzled by the assumption that

[OP] $$ \frac{dk_{2}}{dC_{X}} = 0$$

Why would the steady state concentration be independent of one of the rate constants? That is counter-intuitive and requires some explanation.

Answer 2

[my comment] There are two stable equilibria. There are details and a clear discussion in MIra et al. J. Chem. Educ. 2003, v80, p1488, also a v. thorough explanation in D. Gillespie, 'Markov Processes' publ Academic Press 1992.

You have to solve the cubic equation to find the deterministic values for the equilibrium values, or numerically integrate to get the time profile, but then $k_1$ is divided by 2 and $k_2$ by 6 to relate the stochastic calculation to the usual rate equation. This is because the propensities can be approximated when the number of molecules becomes v. large.

Using the values given in the answer by Karsten the polynomial to solve is

$$\displaystyle ax^3+bx^2+cx+d=0,\quad a=k_2/6, b=-k_1A0/2, c=k_4,d=-k_3B0$$

which has roots at $566.9, 247.6, 85.5$. The middle value is a 'barrier' that the deterministic equations (i.e. normal rate eqns) cannot cross. (When the parameters as such that they produce only one real root then only one equilibrium will exist).

You can see the three roots in the plot below as horizontal lines.

The blue horizontal lines are the roots and the grey lines the rate equation solutions starting at X0 values $0, 50, 150, 240, 260, 450, 700$ which you can see in the plot. You can see how the population is 'attracted' to the upper or lower root but never cross the barrier root at $247$.

The stochastic calculation is very different as there is a chance to rise or fall to an equilibrium value (a root of the cubic) and the barrier can be crossed or recrossed as at each step as there is a random probability of going one way or the other.

The second plot shows this and I have chosen to illustrate that it behaves much as the deterministic plot does, but not all do as can be seen the the third plot.

Having multiple equilibria is of course a result of having a cubic equation. Thus this is unrealistic chemically as a termolecular reaction has such a small chance of occurring that most termolecular reactions are found experimentally to involve two steps.