# What are the possible concentration pathways for a 3-species system?

Inspired by this question, which asks if a multi-species system could oscillate around its equilibrium, I thought about the possible concentration pathways the following more narrow specified system could take:

$$\ce{A ->[\text{k_1}] B}$$ $$\ce{B ->[\text{k_2}] C}$$ $$\ce{C ->[\text{k_3}] A}$$

What I was interested in are the ratio of those concentrations $A:B:C$. The initial conditions are $A(0)=1$, $B(0)=0$, $C(0)=0$. (For simplicity, $[\ce{A}] \equiv A$, and analogous for $B$ and $C$)

What are the possible concentration pathways for different values of $k_1$, $k_2$ and $k_3$?

We have a cyclic cascade of unimolecular chemical reactions

$$\ce{A} \xrightarrow{\kappa_1} \ce{B} \xrightarrow{\kappa_2} \ce{C} \xrightarrow{\kappa_3} \ce{A}$$

where $\kappa_1, \kappa_2, \kappa_3 > 0$ are the rate constants. Assuming that each chemical reaction in this cascade has mass action kinetics, then we have the following system of linear ODEs

$$\begin{bmatrix} \dot a\\ \dot b\\ \dot c\end{bmatrix} = \underbrace{\begin{bmatrix} - \kappa_1 & 0 & \kappa_3\\ \kappa_1 & - \kappa_2 & 0\\ 0 & \kappa_2 & - \kappa_3\end{bmatrix}}_{=: \mathrm M} \begin{bmatrix} a\\ b\\ c\end{bmatrix}$$

where $a := [\ce{A}]$, $b := [\ce{B}]$ and $c := [\ce{C}]$ are the concentrations of the species. Since $\dot a + \dot b + \dot c = 0$, integrating we obtain the following conservation law

$$a (t) + b (t) + c (t) = a_0 + b_0 + c_0$$

where $a_0, b_0, c_0 \geq 0$ are the initial concentrations.

### Stability

The characteristic polynomial of matrix $\mathrm M$ is

$$\det (s \mathrm I_3 - \mathrm M) = s \underbrace{\left( s^2 + (\kappa_1 + \kappa_2 + \kappa_3) s + (\kappa_1 \kappa_2 + \kappa_1 \kappa_3 + \kappa_2 \kappa_3) \right)}_{=: q (s)} = s \, q(s)$$

Since $\kappa_1, \kappa_2, \kappa_3 > 0$, all three coefficients of the (monic) quadratic polynomial $q$ are positive. Hence, both roots of $q$ have negative real parts and, thus, the cyclic cascade is (internally) marginally stable. Given nonzero initial concentrations, the concentrations will not grow unbounded. However, due to the zero eigenvalue, the concentrations will not converge to zero either.

If the concentrations of the species do not converge to zero, then to what values do they converge? At steady state, $\dot a = \dot b = \dot c = 0$. Hence, we have the following homogeneous linear system

$$\begin{bmatrix} - \kappa_1 & 0 & \kappa_3\\ \kappa_1 & - \kappa_2 & 0\\ 0 & \kappa_2 & - \kappa_3\end{bmatrix} \begin{bmatrix} \bar a\\ \bar b\\ \bar c\end{bmatrix} = \begin{bmatrix} 0\\ 0\\ 0\end{bmatrix}$$

where $\bar a, \bar b, \bar c$ are the steady state concentrations. Using Gaussian elimination, we eventually obtain a line passing through the origin that is parameterized as follows

$$\begin{bmatrix} \bar a\\ \bar b\\ \bar c\end{bmatrix} \in \left\{ \gamma \begin{bmatrix} \frac{1}{\kappa_1}\\ \frac{1}{\kappa_2}\\ \frac{1}{\kappa_3}\end{bmatrix} : \gamma \in \mathbb R \right\}$$

where

$$\tau_1 := \frac{1}{\kappa_1} \qquad\qquad\qquad \tau_2 := \frac{1}{\kappa_2} \qquad\qquad\qquad \tau_3 := \frac{1}{\kappa_3}$$

are time constants. Intersecting this line with the plane $\bar a + \bar b + \bar c = a_0 + b_0 + c_0$, we obtain

$$\gamma = \left( \dfrac{1}{\tau_{1} + \tau_{2} + \tau_{3}}\right) \left( a_0 + b_0 +c_0 \right)$$

and, thus, the steady state concentrations are

$$\boxed{\begin{array}{rl} & \\ \qquad\bar a &= \left( \dfrac{\tau_{1}}{\tau_{1} + \tau_{2} + \tau_{3}}\right) \left( a_0 + b_0 +c_0 \right) \qquad \\\\ \qquad\bar b &= \left( \dfrac{\tau_{2}}{\tau_{1} + \tau_{2} + \tau_{3}}\right) \left( a_0 + b_0 +c_0 \right) \qquad\\\\ \qquad\bar c &= \left( \dfrac{\tau_{3}}{\tau_{1} + \tau_{2} + \tau_{3}}\right) \left( a_0 + b_0 +c_0 \right)\\ & \end{array}}$$

If, say, $\kappa_1 \gg \kappa_2, \kappa_3$, i.e., if $\tau_1 \ll \tau_2, \tau_3$, then $\bar a \ll \bar b, \bar c$. If $\kappa_1 = \kappa_2 = \kappa_3$, then $\bar a = \bar b = \bar c$.

Assuming that at least one of $a_0, b_0, c_0$ is nonzero, the ratios of the steady state concentrations are

$$\frac{\bar a}{\bar b} = \frac{\kappa_2}{\kappa_1} = \frac{\tau_1}{\tau_2} \qquad\qquad\qquad \frac{\bar a}{\bar c} = \frac{\kappa_3}{\kappa_1} = \frac{\tau_1}{\tau_3} \qquad\qquad\qquad \frac{\bar b}{\bar c} = \frac{\kappa_3}{\kappa_2} = \frac{\tau_2}{\tau_3}$$

Hence, the ratio of the steady state concentrations of two species is the ratio of their time constants.

### Oscillations

Let

$$\Delta (\kappa_1, \kappa_2, \kappa_3) := (\kappa_1 + \kappa_2 + \kappa_3)^2 - 4 (\kappa_1 \kappa_2 + \kappa_1 \kappa_3 + \kappa_2 \kappa_3) \\ = (\kappa_1 - \kappa_2)^2 + (\kappa_1 - \kappa_3)^2 + (\kappa_2 - \kappa_3)^2 - (\kappa_1^2 + \kappa_2^2 + \kappa_3^2)$$

be the discriminant of quadratic polynomial $q$. When $\Delta < 0$, the roots of $q$ have nonzero imaginary parts and, thus, the cyclic cascade will exhibit damped oscillations. Hence, when

$$\kappa_1^2 + \kappa_2^2 + \kappa_3^2 > (\kappa_1 - \kappa_2)^2 + (\kappa_1 - \kappa_3)^2 + (\kappa_2 - \kappa_3)^2$$

the cyclic cascade will exhibit oscillatory behavior. Provided that the inequality above is satisfied, the concentrations will oscillate as follows

$$\begin{array}{rl} a (t) &= \bar a + \eta_A \, e^{- \kappa_0 t} \cos( \omega_0 t - \phi_A)\\ b (t) &= \bar b + \eta_B \, e^{- \kappa_0 t} \cos( \omega_0 t - \phi_B)\\ c (t) &= \bar c + \eta_C \, e^{- \kappa_0 t} \cos( \omega_0 t - \phi_C)\end{array}$$

where

$$\begin{array}{rl} \kappa_0 &:= \frac{1}{2} \left( \kappa_1 + \kappa_2 + \kappa_3 \right)\\ \omega_0 &:= \frac 12 \sqrt{\kappa_1^2 + \kappa_2^2 + \kappa_3^2 - (\kappa_1 - \kappa_2)^2 - (\kappa_1 - \kappa_3)^2 - (\kappa_2 - \kappa_3)^2}\end{array}$$

Each species will have its own magnitude and phase lag that are dependent on the initial conditions.

• Nice to see a more complete answer. – John H. K. Oct 18 '17 at 5:50
• @JohnH.K. I updated my answer. I think the most interesting part of it is the fact that the steady state concentrations of $\ce{A}$, $\ce{B}$ and $\ce{C}$ are proportional to $\frac{1}{\kappa_1}$, $\frac{1}{\kappa_2}$ and $\frac{1}{\kappa_3}$, respectively. This agrees with your diagrams. – Rodrigo de Azevedo Oct 18 '17 at 18:20

Of course, we can't display every pathway there is. We have to select a wide variety of rate constants to get a "good map" of the parameter-space. The approach I took was this:

1. Construct the differential equations

$$\frac{\mathrm dA}{\mathrm dt} = -k_1 A + k_3 C$$ $$\frac{\mathrm dB}{\mathrm dt} = -k_2 B + k_1 A$$ $$\frac{\mathrm dC}{\mathrm dt} = -k_3 C + k_2 B$$

1. Solve them numerically for every possible combination of 0.1, 1 and 10 for each rate constant and the above stated initial conditions
2. This leads to 27 pathways, which are plotted in a ternary diagram

The title of the ternary diagram is $k_1:k_2:k_3$, which is 0.1 : 0.1 : 0.1 in the first one.

(click for larger image)

• What software did you use to generate these diagrams? – Rodrigo de Azevedo Oct 19 '17 at 7:35
• @RodrigodeAzevedo Honestly Im not completely sure, either Python / Matplotlib or R. Maybe I'll find the source code when I'm home. – John H. K. Oct 19 '17 at 8:43