Test Models for Statistical Inference: Two-Dimensional Reaction Systems Displaying Limit Cycle Bifurcations and Bistability

Abstract

Theoretical results regarding two-dimensional ordinary-differential equations (ODEs) with second-degree polynomial right-hand sides are summarized, with an emphasis on limit cycles, limit cycle bifurcations, and multistability. The results are then used for construction of two reaction systems, which are at the deterministic level described by two-dimensional third-degree kinetic ODEs. The first system displays a homoclinic bifurcation, and a coexistence of a stable critical point and a stable limit cycle in the phase plane. The second system displays a multiple limit cycle bifurcation, and a coexistence of two stable limit cycles. The deterministic solutions (obtained by solving the kinetic ODEs) and stochastic solutions [noisy time-series generating by the Gillespie algorithm, and the underlying probability distributions obtained by solving the chemical master equation (CME)] of the constructed systems are compared, and the observed differences highlighted. The constructed systems are proposed as test problems for statistical methods, which are designed to detect and classify properties of given noisy time-series arising from biological applications.

Appendices

Appendix 1: Perturbed x-Factorable Transformation

Definition 1

Consider applying an x-factorable transformation, as defined in [2], on (5), and then adding to the resulting right-hand side a zero-degree term ɛ v, with ɛ ≥ 0 and vector v = (1, 1)^⊤, resulting in

$$\displaystyle\begin{array}{rcl} \frac{\mathrm{d}\mathbf{x}} {\mathrm{d}t} & =\varepsilon \mathbf{v} + \mathcal{X}(\mathbf{x})\boldsymbol{\mathcal{P}}(\mathbf{x};\,\mathbf{k})\, =\,\varepsilon \mathbf{v} + (\Psi _{\mathcal{X}}\boldsymbol{\mathcal{P}})(\mathbf{x};\,\mathbf{k}) \equiv (\Psi _{\mathcal{X}_{\varepsilon }}\boldsymbol{\mathcal{P}})(\mathbf{x};\,\mathbf{k}).&{}\end{array}$$

(20)

Then $\Psi _{\mathcal{X}_{\varepsilon }}: \mathbb{P}_{2}(\mathbb{R}^{2};\, \mathbb{R}^{2}) \rightarrow \mathbb{P}_{3}(\mathbb{R}^{2};\, \mathbb{R}^{2})$, mapping $\boldsymbol{\mathcal{P}}(\mathbf{x};\,\mathbf{k})$ to $(\Psi _{\mathcal{X}_{\varepsilon }}\boldsymbol{\mathcal{P}})(\mathbf{x};\,\mathbf{k})$, is called a perturbed x-factorable transformation if ɛ ≠ 0. If ɛ = 0, the transformation reduces to an (unperturbed) x-factorable transformation, $\Psi _{\mathcal{X}}\equiv \Psi _{\mathcal{X}_{0}}$, defined in [2].

Lemma 1

$(\Psi _{\mathcal{X}_{\varepsilon }}\boldsymbol{\mathcal{P}})(\mathbf{x};\,\mathbf{k})$ from Definition 1 is a kinetic function, i.e. $(\Psi _{\mathcal{X}_{\varepsilon }}\boldsymbol{\mathcal{P}})(\mathbf{x};\,\mathbf{k}) \in \mathbb{P}_{3}^{\mathcal{K}}(\mathbb{R}_{\geq }^{2};\, \mathbb{R}^{2})$ .

Proof

$(\Psi _{\mathcal{X}}\boldsymbol{\mathcal{P}})(\mathbf{x};\,\mathbf{k})$ is a kinetic function [2]. Since, from (20), $(\Psi _{\mathcal{X}_{\varepsilon }}\boldsymbol{\mathcal{P}})(\mathbf{x};\,\mathbf{k}) =\varepsilon \mathbf{v} + (\Psi _{\mathcal{X}}\boldsymbol{\mathcal{P}})(\mathbf{x};\,\mathbf{k})$, with ɛ ≥ 0 and v = (1, 1)^⊤, it follows that $(\Psi _{\mathcal{X}_{\varepsilon }}\boldsymbol{\mathcal{P}})(\mathbf{x};\,\mathbf{k})$ is kinetic as well. □

We now provide a theorem relating location, stability and type of the positive critical points of (5) and (20).

Theorem 1

Consider the ODE system (5) with positive critical points $\mathbf{x}^{{\ast}} \in \mathbb{R}_{>}^{2}$ . Let us assume that $\mathbf{x}^{{\ast}}\in \mathbb{R}_{>}^{2}$ is hyperbolic, and is not the degenerate case between a node and a focus, i.e. it satisfies the condition

$$\displaystyle\begin{array}{rcl} \left (\mathit{\text{tr}}\left (\nabla \boldsymbol{\mathcal{P}}(\mathbf{x}^{{\ast}};\,\mathbf{k})\right )\right )^{2} - 4\mathit{\text{det}}\left (\nabla \boldsymbol{\mathcal{P}}(\mathbf{x}^{{\ast}};\,\mathbf{k})\right )& \neq 0,&{}\end{array}$$

(21)

as well as conditions (ii) and (iii) of Theorem 3.3 in [2]. Then positivity, stability and type of the critical point $\mathbf{x}^{{\ast}} \in \mathbb{R}_{>}^{2}$ are invariant under the perturbed x-factorable transformations $\Psi _{\mathcal{X}_{\varepsilon }}$ , for sufficiently small ɛ ≥ 0. Assume (5) does not have boundary critical points. Consider the two-dimensional ODE system (20) with ɛ = 0, and with boundary critical points denoted $\bar{\mathbf{x}}^{0} \in \mathbb{R}_{\geq }^{2}$ , $\bar{\mathbf{x}}^{0} = (\bar{x}_{b,1}^{0},\bar{x}_{b,2}^{0})$ , $\bar{x}_{b,1}^{0}\bar{x}_{b,2}^{0} = 0$ . Assume that for i ∈ {1, 2}

$$\displaystyle\begin{array}{rcl} \frac{\partial \mathcal{P}_{i}(\mathbf{\bar{x}}_{b}^{0};\,\mathbf{k})} {\partial x_{i}} & \neq 0,\,\,\,\,\,\,\,\,\,\mathit{\text{if }}\,\,\,\bar{x}_{b,i}^{0}\,\neq \,0,&{}\end{array}$$

(22)

and that for some i ∈ {1, 2}

$$\displaystyle\begin{array}{rcl} \mathcal{P}_{i}(\mathbf{\bar{x}}_{b}^{0};\,\mathbf{k})&> 0,\,\,\,\,\,\,\,\,\,\mathit{\text{if }}\,\,\,\bar{x}_{ b,i}^{0}\, =\, 0.&{}\end{array}$$

(23)

Then, the critical point $\mathbf{\bar{x}}_{b}^{0} \in \mathbb{R}_{\geq }^{2}$ of the two-dimensional ODE system (20) with ɛ = 0 becomes the critical point $\mathbf{\bar{x}}_{b}\notin \mathbb{R}_{\geq }^{2}$ of system (20) for sufficiently small ɛ > 0.

Proof

The critical points of (20) are solutions of the following regularly perturbed algebraic equation

$$\displaystyle\begin{array}{rcl} \varepsilon \mathbf{v} + \mathcal{X}(\mathbf{\bar{x}})\boldsymbol{\mathcal{P}}(\mathbf{\bar{x}};\,\mathbf{k})& = \mathbf{0}.&{}\end{array}$$

(24)

Let us assume $\mathbf{\bar{x}}$ can be written as the power series

$$\displaystyle\begin{array}{rcl} \mathbf{\bar{x}}& = \mathbf{\bar{x}}^{0} +\varepsilon \mathbf{\bar{x}}^{1} + \mathcal{O}(\varepsilon ^{2}),&{}\end{array}$$

(25)

where $\mathbf{\bar{x}}^{0} \in \mathbb{R}_{\geq }^{2}$ are the critical points of (20) with ɛ = 0. Substituting the power series (25) into (24), and using the Taylor series theorem on $\boldsymbol{\mathcal{P}}(\mathbf{\bar{x}};\,\mathbf{k})$, so that $\boldsymbol{\mathcal{P}}(\mathbf{\bar{x}}^{0} +\varepsilon \mathbf{\bar{x}}^{1} + \mathcal{O}(\varepsilon ^{2});\,\mathbf{k}) =\boldsymbol{ \mathcal{P}}(\mathbf{\bar{x}}^{0};\,\mathbf{k}) +\varepsilon \nabla \boldsymbol{\mathcal{P}}(\mathbf{\bar{x}}^{0};\,\mathbf{k})\mathbf{\bar{x}}^{1} + \mathcal{O}(\varepsilon ^{2})$, as well as that $\mathcal{X}(\mathbf{\bar{x}}) = \mathcal{X}(\mathbf{\bar{x}}^{0}) +\varepsilon \mathcal{X}(\mathbf{\bar{x}}^{1}) + \mathcal{O}(\varepsilon ^{2})$, and equating terms of equal powers in ɛ, the following system of polynomial equations is obtained:

$$\displaystyle\begin{array}{rcl} & & \mathcal{O}\left (1\right ):\; \mathcal{X}(\mathbf{\bar{x}}^{0})\boldsymbol{\mathcal{P}}(\mathbf{\bar{x}}^{0};\,\mathbf{k})\, = 0, \\ & & \mathcal{O}(\varepsilon ):\; \mathcal{X}(\mathbf{\bar{x}}^{0})\nabla \boldsymbol{\mathcal{P}}(\mathbf{\bar{x}}^{0};\,\mathbf{k})\mathbf{\bar{x}}^{1} + \mathcal{X}(\mathbf{\bar{x}}^{1})\boldsymbol{\mathcal{P}}(\mathbf{\bar{x}}^{0};\,\mathbf{k})\, = -\mathbf{v}.{}\end{array}$$

(26)

Order 1 equation. The positive critical points $\mathbf{\bar{x}}^{0} \in \mathbb{R}_{>}^{2}$ satisfy $\boldsymbol{\mathcal{P}}(\mathbf{\bar{x}}^{0};\,\mathbf{k}) = \mathbf{0}$. Since $\boldsymbol{\mathcal{P}}(\mathbf{x};\,\mathbf{k})$ has no boundary critical points by assumption, critical points $\mathbf{\bar{x}}_{b}^{0} \in \mathbb{R}_{\geq }^{2}$ with $\bar{x}_{b,i}^{0} = 0$, $\bar{x}_{b,j}^{0}\neq 0$, $\bar{x}_{b,1}^{0}\bar{x}_{b,2}^{0} = 0$, i, j ∈ {1, 2}, satisfy $\mathcal{P}_{i}(\mathbf{\bar{x}}_{b}^{0};\,\mathbf{k})\neq 0$, $\mathcal{P}_{j}(\mathbf{\bar{x}}_{b}^{0};\,\mathbf{k}) = 0$.

Order ɛ equation. Vector $\mathbf{\bar{x}}^{1}$, corresponding to a positive $\mathbf{\bar{x}}^{0}$, satisfies

$$\displaystyle\begin{array}{rcl} \mathcal{X}(\mathbf{\bar{x}}^{0})\nabla \boldsymbol{\mathcal{P}}(\mathbf{\bar{x}}^{0};\,\mathbf{k})\mathbf{\bar{x}}^{1}& = -\mathbf{v},& {}\\ \end{array}$$

which can be solved provided $\mathbf{\bar{x}}^{0}$ is a hyperbolic critical point. Vector $\mathbf{\bar{x}}_{b}^{1}$, corresponding to a nonnegative $\mathbf{\bar{x}}_{b}^{0}$, is given by

$$\displaystyle\begin{array}{rcl} \bar{x}_{b,i}^{1} = \left \{\begin{array}{@{}l@{\quad }l@{}} -(\mathcal{P}_{i}(\mathbf{\bar{x}}_{b}^{0};\,\mathbf{k}))^{-1}, \quad &\text{if }\bar{x}_{b,i}^{0} = 0, \\ \left (\frac{\partial \mathcal{P}_{i}(\mathbf{\bar{x}}_{b}^{0};\,\mathbf{k})} {\partial x_{i}} \right )^{-1}\left ((\mathcal{P}_{j}(\mathbf{\bar{x}}_{b}^{0};\,\mathbf{k}))^{-1}\frac{\partial \mathcal{P}_{i}(\mathbf{\bar{x}}_{b}^{0};\,\mathbf{k})} {\partial x_{j}} - (\bar{x}_{b,i}^{0})^{-1}\right ),\quad &\text{if }\bar{x}_{b,i}^{0}\neq 0, \end{array} \right.& & {}\\ \end{array}$$

from which conditions (22) and (23) follow. □

Appendix 2: Bicyclic System with Large Attractors

Consider the following deterministic kinetic equations

$$\displaystyle\begin{array}{rcl} \frac{\mathrm{d}x_{1}} {\mathrm{d}t} & =& k_{1} + x_{1}(-k_{2} + k_{3}x_{1} + k_{4}x_{2} - k_{5}x_{1}x_{2}), \\ \frac{\mathrm{d}x_{2}} {\mathrm{d}t} & =& k_{6} + x_{2}(k_{7} - k_{8}x_{1} + k_{9}x_{2} + k_{10}x_{1}^{2} - k_{ 11}x_{2}^{2}),{}\end{array}$$

(27)

with the coefficients k given by

$$\displaystyle\begin{array}{rcl} k_{1}& =& 10^{-3},\;\;\;k_{ 2} = 10,\;\;\;k_{3} = 1,\;\;\;k_{4} = 1,\;\;\;k_{5} = 0.1,\;\;\;k_{6} = 10^{-3}, \\ k_{7}& =& 3.7,\;\;\;k_{8} = 1.9,\;\;\;k_{9} = 1.01,\;\;\;k_{10} = 0.1,\;\;\;k_{11} = 0.05. {}\end{array}$$

(28)

The canonical reaction network induced by system (27), involving two species s ₁ and s ₂ and eleven reactions r ₁, r ₂, …, r ₁₁ under mass-action kinetics, is given by

(29)

In Fig. 5a, we show the two stable limit cycles obtained by numerically solving (27) with parameters (28). In Fig. 5b, in addition to the limit cycles, we also show in blue a representative sample path obtained by applying the Gillespie algorithm on (29). Let us note that (27) was constructed in a similar fashion as system (12) in Sect. 3.2, using the results from [40, 46].

Fig. 5

Panel (a) displays numerically approximated stable limit cycles L ₁ and L ₃ in the state-space of system (27), with parameters (28) and reactor volume V = 100. Panel (b) displays in blue a representative sample path, generated by applying the Gillespie algorithm on the underlying reaction network (29) for the same parameters as in panel (a). Also shown are two deterministic trajectories, one initiated near the limit cycle L ₁, while the other near L ₃. One can observe that the stochastic sample path switches between the two deterministic attractors