# Asymmetric autocatalytic reactions and their stationary distribution

Cameron Gallinger, Lea Popovic Department of Mathematics and Statistics
Concordia University, Montreal QC H3G 1M8, Canada; email:
<EMAIL_ADDRESS>
(November 27. 2023)

###### Abstract

We consider a general class of autocatalytic reactions, that has been shown to
display stochastically switching behaviour (Discreteness Induced Transitions)
in some parameter regimes. This behaviour was shown to occur when either the
overall species count is low, or when the rate of inflow and outflow of
species is relatively much smaller than the rate of autocatalytic reactions.
The long-term behaviour of this class was analyzed in [20] with an analytic
formula for the stationary distribution in the symmetric case. We focus on the
case of asymmetric autocatalytic reactions and provide a formula for an
approximate stationary distribution of the model. We show this distribution
has different properties corresponding to the distinct behaviour of the
process in the three parameter regimes. In the Discreteness Induced
Transitions regime, the formula provides the fraction of time spent at each of
the stable points.

†† _Keywords and phrases._ reaction networks, autocatalytic reactions,
stationary distribution, discreteness induced transitions, Moran model, genic
selection

##  1 Introduction

Many important processes in biology rely on the switching behaviour between
distinct states of the internal state of the system. Mathematically, this may
be the result of alternating between different stable states for the dynamics
of the molecular composition of the system. This feature is present in many
gene-expression systems, where a gene alternates between two types of states
(“on” and “off”) regulating the production of a protein. It is also present in
many phosphorylation switches in signalling pathways. Such bistable switching
patterns arise from distinct stable equilibria in the deterministic dynamics
and the ability of infrequent large stochastic fluctuations to pull the system
from a basin of attraction of one equilibrium to the other. However, there are
also cases of chemical dynamics in which bistability is not possible in the
deterministic model and is only possible in the stochastic model of the same
chemical reaction system (e.g. [5]). The importance of stochastic effects in
biological switching is well recognized and a variety of stochastic processes
are being used to model them (see [1] for a review).

Autocatalytic reaction networks are an important class of reaction networks
whose long-term behaviour can substantively differ when modelled
deterministically versus stochastically (e.g. [6]). Autocatalytic reactions
are of broad interest in modelling living systems as they involve only a small
number of distinct species and interactions, but they exhibit a range of
outcomes such as self-sustaining growth, oscillations and symmetry breaking.
Autocatalysis occurs in many elements of cellular metabolism including
glycolysis, mitosis, apoptosis, and DNA replication (for a review see [3]),
and may plausibly play a key role in the origin of life (e.g.[4]). They were
recently used in chemical computing, [2], as part of an artificial neural
network design to perform digital image recognition tasks.

We analyze a model for a simple example of an autocatalytic reaction network
proposed by Togashi-Kaneko in [18]. They were particularly interested in the
dynamics of the network when the total species count in the system was low,
and observed numerically that stochasticity at the low species count led to a
new type of switching dynamics. The system spent most of its time in states
where some combination of species is absent, switching rapidly by a successive
increase in some species count from one such state to another. They called
this behaviour Discretely Induced Transitions (DITs) because discrete changes
in species counts drive the switching. This dynamics is prominent in the
parameter regime when the rates of autocatalytic reactions are fast relative
to the rates of inflow and outflow from a reservoir, and the average time
between switches is large.

The long-term behaviour of this model was further explored in [7], where the
same authors observed that, if only the long term averages are observed,
autocatalytic reaction networks with strong symmetries do not show the
important distinction between small and large species counts. However, in
systems possessing asymmetry, numerical simulations clearly show that the
shape of the long-term distributions shifts as the species count decreases.
They observed similar effects in simple catalytic reaction networks in [9],
where the low count of some species type effectively acts by switching on and
off parts of the reaction network. They also observed this effect in spatial
models in [8], where higher abundances are typically modelled by reaction-
diffusion systems, but the localization and low abundance of some species
types generates different possibilities for the long term behaviour patterns
of the system.

First analytic results for symmetric auto-catalytic reaction networks were
obtained by a stochastic differential equation approximation of the system
dynamics in [10]. They used a scaling parameter composed of the product of the
system volume and inflow rate, and derived an analytic expression for the
stationary distribution of a linear combination of the species in the system,
based on a time-scale separation. In [11] they used the same tools to explore
a version of this auto-catalytic reaction network without inflows and
outflows, but with single species conversions. They provided a formula for the
mean switching time between the states when the system contains only one of
the species. A rigorous approximation by an obliquely reflected stochastic
differential equation of a general autocatalytic reaction network was
subsequently developed in [12], where a constrained Langevin approximation
(based on the theory of Leite and Williams [13]) was analyzed, together with
the stationary distribution of this reflected stochastic differential
equation.

Stochastic approximations are sensitive to the different abundances in
molecular amounts (e.g. [14]), as well as to the time-scale separations in the
rates of different reactions in the network (e.g. [15]). For this reason, it
is worth analyzing the properties of reaction networks from the Markov chain
model directly. Analytic results for the Markov chain model were explored in
[19], where the effect of species counts was shown on the flow of the
molecular ’current’ from one species type to another. Analytic results for the
stationary distribution of a related class of models were derived in [16], for
reaction networks without inflows and outflows and hence conserved overall
species counts. A thorough analytic exploration of the auto-catalytic reaction
network with inflows and outflows was performed in [20], where it was shown
that the Markov chain model of a general auto-catalytic reaction network is
exponentially ergodic. The authors further analyzed the symmetric version of
this model and derived an explicit form for the stationary distribution as a
Dirichlet-multinomial distribution (Theorem 4.3, [20]). However, an explicit
formula for the case of asymmetric autocatalytic reaction networks remained
unknown.

To analyze the stationary distribution in the asymmetric case, we explored the
connection of the auto-catalytic model with inflows and outflows to that of
the Moran model with mutation from population genetics. A connection with the
Moran model was noted in [11, 19] for symmetric rates of auto-catalysis. We
observe that asymmetric rates appear in Moran models that contain genic
selection, as in [21], and the formula for the stationary distribution for
such a Moran model is the starting point in our analysis. We explore the
consequences due to the differences between the two models and explore how
close the stationary distribution for the Moran model is to that of the auto-
catalytic network. We analyze the effect of the asymmetry in the long-term
distribution of the species count in the system. We show the approximate
stationary distribution also has the signature of the Discretely Induced
Transitions in the parameter regime corresponding to low species counts, or
slow inflow and outflow rates in the system.

##  2 Stochastic (Markov chain) model for auto-catalytic reactions

The auto-catalytic reaction model we examine is the Togashi-Kaneko model for
an open system. It consists of $d$ molecular species
$\\{A_{i}\\}_{i=1,\dots,d}$ involved in auto-catalytic reactions and subject
to inflows and outflows from the system:

$A_{i}+A_{j}\xrightarrow{\kappa_{i}}2A_{i}\quad\quad\emptyset\xrightleftharpoons[\delta]{\lambda_{i}}A_{i}\qquad
i,j=1,2,\dots,d$ (2.1)

The stochastic model is a continuous-time Markov chain $(X(t))_{t\geq 0}$,
where $X(t)=(X_{1}(t),\\!\dots\\!,X_{d}(t))$ counts the number of molecules of
species $i\in\\{1,\dots,d\\}$ present in the system at time $t\geq 0$. The
state space for the Markov chain is
$E=\\{a=(a_{1},\dots,a_{d})\in\mathbb{N}^{d}\\}$ and the transition rates are
given by mass-action kinetics as

$\begin{split}&q_{a,a-e_{j}+e_{i}}=\kappa_{i}a_{i}a_{j},\qquad\forall i\neq
j\in\\{1,\dots,d\\}\\\ &q_{a,a+e_{i}}=\lambda\quad q_{a,a-e_{i}}=\delta
a_{i},\quad\forall i\in\\{1,\dots,d\\}\end{split}$ (2.2)

Note that this is not necessarily a cyclic model (as all pairs of different
species can catalyze each other), and the rates of autocatalytic reactions are
different for each pair. The only constraint on the reactions that we impose
is that the outflow rate is the same for all species, which ensures the
lumpability of the Markov chain into a process that counts the total mass in
the system.

The concept of lumpability means that the state space $E$ can be partitioned
into subsets $E_{n}=\\{a\in\mathbb{N}^{d}|\sum_{i=1}^{d}=n\\}$, and that the
lumped process $(N(t))_{t\geq 0}$ defined by $N(t)=\sum_{i=1}^{d}X_{i}(t)$ is
a continuous time Markov chain on $E_{n}$ with associated transition rates

$q_{n,n+1}=\sum_{i=1}^{d}\lambda_{i}\quad\quad q_{n,n-1}=n\delta.$ (2.3)

Because the total mass process $(N(t))_{t\geq 0}$ is a linear death process
with immigration on $\mathbb{Z}_{+}$, we know that it is irreducible and
positive recurrent with stationary distribution $\nu$ that is Poisson with
mean equal to $\sum_{i=1}^{d}\lambda_{i}/\delta$.

In [20] Theorem 4.1, Bibbona, Kim and Wiuf showed that the Markov chain
$(X(t))_{t\geq 0}$ corresponding to (2.2) is positive recurrent on $E$ with a
unique stationary distribution $\Pi$, and that it converges exponentially fast
to this distribution.

In the special case when all the auto-catalytic rates are the same:
$\kappa_{i}=\kappa,\forall i$, they derive an explicit formula for the
stationary distribution $\Pi$ in the form

$\Pi(a)=\pi(a|n)\nu(n)$ (2.4)

where $\nu\sim$ Poisson$(\mu)$ is the stationary distribution of the total
mass process, and $\pi(a|n)$ is the Dirichlet-multinomial$(n,\alpha)$
distribution on $E_{n}$

$\pi(a|n)={n\choose
a}\frac{\Gamma\big{(}\sum_{i=1}^{d}\alpha_{i}\big{)}}{\Gamma\big{(}n+\sum_{i=1}^{d}\alpha_{i}\big{)}}\prod_{i=1}^{d}\frac{\Gamma\big{(}a_{i}+\alpha_{i}\big{)}}{\Gamma\big{(}\alpha_{i}\big{)}}$
(2.5)

with parameters $\alpha=(\alpha_{1},\dots,\alpha_{d})$ given by

$\alpha_{i}=\frac{\delta\lambda_{i}}{\kappa\sum_{j=1}^{d}\lambda_{j}}$ (2.6)

This analytic formula for $\Pi$ from (2.4) and (2.5) is directly verified
using the time evolution for the probability distribution of the system.

###  2.1 Connection to population genetics

The Moran model is a fundamental stochastic model for the evolution of a
population of species types undergoing reproduction and mutation in continuous
time (overlapping generations). The model describes a population of fixed size
$n$, with $d$ different species types $\\{A_{i}\\}_{i=1,\dots,d}$ , whose
number of individuals at time $t\geq 0$ is denoted by
$\tilde{X}(t)=(\tilde{X}_{1}(t),\dots,\tilde{X}_{d}(t))$. Each individual of
type $i$ gives birth at the rate $\kappa_{i}$, and its offspring replaces an
individual chosen at random to die. The species undergo mutations with each
individual of type $i$ changing to type $j$ at rate of $vp_{ij}$, where
$p_{ij}$ are the entries of a stochastic matrix $P$ and $v>0$ is the overall
mutation rate.

The species reproduction and mutation changes can be represented graphically
by:

$A_{i}+A_{j}\xrightarrow{\kappa_{j}}2A_{j}\quad\quad
A_{i}\xrightarrow{vp_{j}}A_{j}$ (2.7)

where we assumed parent-independent mutations, $p_{ij}=p_{j}$. The fact that
the fitness rates $\kappa_{j}$ are not equal across species is due to genic
selection, and the special case $\kappa_{i}=\kappa,\forall
i\in\\{1,\dots,d\\}$ corresponds to a model that is neutral under selection.

The continuous time Markov chain of species counts $(\tilde{X}(t))_{t\geq 0}$
has transition rates from $a$ to $a-e_{i}+e_{j}$ equal to

$\tilde{q}_{a,a-e_{i}+e_{j}}=a_{i}\bigg{[}\frac{\kappa_{j}}{n}a_{j}+\nu
p_{j}\bigg{]}.$ (2.8)

The Moran model differs from the model of auto-catalytic reactions (2.1) in
that it has a conservation of mass (closed system) and lacks the (open system)
effects of the inflow and outflow reactions.

A comparison of long-term behaviour for closed versus open systems exists for
deterministic models of reaction networks: if the open system is created from
the closed system by adding inflow and outflow for each species type, then the
open system has the same number of steady states and stability as the closed
system ([22, 23, 24]). For stochastic systems, the comparison is different as
the closed system may have absorbing states when the open system does not
affect the existence and the number of modes in the stationary distribution.
The reason the Moran model is remarkably close to the TK model in its long-
term behaviour is due to the fact that one can create a correspondence between
the mutation events in the former and the in-/out-flows in the latter (see
explanation above (3.5)).

##  3 Approximate stationary distribution

The Moran model with genic selection on a population of size $n$ was shown by
Etheridge and Griffiths [21] to have a stationary distribution proportional to

$\pi_{n}(a)\propto\kappa_{1}^{a_{1}}\cdots\kappa_{d}^{a_{d}}{n\choose
a}\frac{\alpha_{1_{(a_{1})}}\cdots\alpha_{d_{(a_{d})}}}{|\alpha|_{(n)}},\;a\in
E_{n}$ (3.1)

where the parameters are given by

$\alpha_{i}=\frac{nvp_{i}}{\kappa_{i}},$ (3.2)

$|\alpha|=\sum_{i=1}^{d}\alpha_{i}$, and $\alpha_{(a)}$ is the ascending
factorial notation (Pochammer function)

$\alpha_{(a)}=\alpha(\alpha+1)\textperiodcentered\textperiodcentered\textperiodcentered(\alpha+a-1).$

The normalizing constant for $\pi_{n}(a)$ is given as the partition function

$u(\alpha,\kappa,n)=\mathbb{E}\bigg{[}\bigg{(}\sum_{i=1}^{d}\kappa_{i}\xi_{i}\bigg{)}^{n}\bigg{]}$
(3.3)

where $\xi=(\xi_{1},\xi_{2},...,\xi_{d})$ has a Dirichlet-
multinomial$(n,\alpha)$ distribution. Writing the ascending factorial in terms
of Gamma functions $\alpha_{(a)}=\Gamma(\alpha+a)/\Gamma(\alpha)$, the
stationary distribution can also be written as

$\pi_{n}(a)\propto\prod_{i=1}^{d}\kappa_{i}^{a_{i}}{n\choose
a}\frac{\Gamma(\sum_{i=1}^{d}\alpha_{i})}{\Gamma(\sum_{i=1}^{d}\alpha_{i}+n)}\prod_{i=1}^{d}\frac{\Gamma(\alpha_{i}+a_{i})}{\Gamma(\alpha_{i})}$
(3.4)

In case of the neutral model, when $\kappa_{i}=\kappa,\forall i=1,\dots,d$,
the partition function becomes
$u(\alpha,\kappa,n)=\kappa^{n}\mathbb{E}\left[\left(\sum_{i=1}^{d}\xi_{i}\right)^{n}\right]=\kappa^{n}$
since $\sum_{i=1}^{d}\xi_{i}=1$. The stationary distribution then reduces to
the Dirichlet-multinomial distribution

$\pi_{n}(a)={n\choose
a}\frac{\Gamma(\sum_{i=1}^{d}\alpha_{i})}{\Gamma(\sum_{i=1}^{d}\alpha_{i}+n)}\prod_{i=1}^{d}\frac{\Gamma(\alpha_{i}+a_{i})}{\Gamma(\alpha_{i})}$

which is the same as the stationary distribution $\pi(a|n)$ in (2.5) derived
by Bibbona, Kim and Wiuf under the parameter change

$\alpha_{i}=\frac{\delta\lambda_{i}}{\kappa\sum_{j=1}^{d}\lambda_{j}}.$

To create a correspondence between the Moran model with genic selection and
the auto-catalytic reactions model, suppose we require that each outflow of
some species is identified in time with the next inflow of some species (e.g.
$A_{i}$), thus closing the system. With this modification the overall rate of
inflow of species $A_{i}$ (simultaneously with some outflow) would be
$\delta\lambda_{i}/(\sum_{j=1}^{d}\lambda_{j})$. In the Moran model the
overall rate of mutation into species $A_{i}$ is $nvp_{i}$. Identifying these
inflow rates in the equation for the parameter (3.2) would result in using the
parameter

$\alpha_{i}=\frac{\delta\lambda_{i}}{\kappa_{i}\sum_{j=1}^{d}\lambda_{j}}$
(3.5)

for the asymmetric auto-catalytic reaction system.

The correspondence suggests the distribution $\pi_{n}$ from (3.4) may be a
good candidate in forming the approximate stationary distribution for
asymmetric autocatalytic reaction system, and we use $\tilde{\pi}_{n}$ to
denote $\pi_{n}$ with parameters $\alpha_{i}$ satisfying (3.5). Our proposed
approximate stationary distribution for the autocatalytic reaction system
(2.1) then is

$\tilde{\Pi}(a)=\nu(n)\frac{1}{u(\alpha,\kappa,n)}\bigg{[}\prod_{i=1}^{d}\kappa_{i}^{a_{i}}\bigg{]}{n\choose
a}\frac{\Gamma(\sum_{i=1}^{d}\alpha_{i})}{\Gamma(\sum_{i=1}^{d}\alpha_{i}+n)}\prod_{i=1}^{d}\frac{\Gamma(\alpha_{i}+a_{i})}{\Gamma(\alpha_{i})},\;\;n=\sum_{i=1}^{d}a_{i}$
(3.6)

We next analyze how well the distribution $\tilde{\Pi}$ approximates the true
stationary distribution $\Pi$.

###  3.1 Analytical results

To assess how close $\tilde{\Pi}$ is to the true stationary distribution of
(2.1), let $\mathcal{A}$ denote the generator of the Markov chain
$(X(t))_{t\geq 0}$ defined for any $f:E\mapsto\mathbb{R}$ by

$\displaystyle{\mathcal{A}}f(a)$
$\displaystyle=\sum_{i,j=1}^{d}q_{a,a-e_{j}+e_{i}}[f(a-e_{j}+e_{i})-f(a)]$
(3.7)
$\displaystyle+\sum_{i=1}^{d}q_{a,a+e_{i}}[f(a+e_{i})-f(a)]+\sum_{i=1}^{d}q_{a,a-e_{i}}[f(a-e_{i})-f(a)]$

where the transition rates $q_{i,j}$ are as in (2.2). Let $p_{t}(\mu,a)$
denote the probability that the process $(X(t))_{t\geq 0}$ is in state $a\in
E$ at time $t\geq 0$, given that its initial distribution is $X(0)\sim\mu$.
Let ${\mathcal{A}}^{*}$ be the adjoint of ${\mathcal{A}}$ defined by

$\displaystyle{\mathcal{A}}^{*}p_{t}(\mu,a)$ $\displaystyle=\sum_{i\neq
j=1}^{d}\big{[}q_{a+e_{j}-e_{i},a}\,p_{t}(\mu,a+e_{j}-e_{i})\mathbf{1}_{a+e_{j}-e_{i}\in
E}-q_{a,a+e_{j}-e_{i}}p_{t}(\mu,a)\big{]}$ (3.8)
$\displaystyle+\sum_{i=1}^{d}\big{[}q_{a+e_{i},a}p_{t}(\mu,a+e_{i})-q_{a,a+e_{i}}p_{t}(\mu,a)\big{]}+\sum_{i=1}^{d}\big{[}q_{a-e_{i},a}p_{t}(\mu,a-e_{i})-q_{a,a-e_{i}}p_{t}(\mu,a)\big{]}$

which specifies the evolution of $p_{t}(\mu,a)$ via the Kolmogorov forward
equation (chemical master equation in the biology literature) as

$\frac{d}{dt}p_{t}(\mu,a)={\mathcal{A}}^{*}p_{t}(\mu,a)\mathbf{1}_{a\in E}.$
(3.9)

The stationary distribution $\Pi$ of (2.1) is in the null space of
${\mathcal{A}}^{*}$, satisfying the global balance equation:
${\mathcal{A}}^{*}\Pi(a)=0$. Using the proposed distribution
$\tilde{\Pi}=\nu(n)\tilde{\pi}_{n}(a)$ equation (3.8) becomes

$\displaystyle{\mathcal{A}}^{*}\tilde{\Pi}(a)$
$\displaystyle=\nu(n)\sum_{i\neq
j=1}^{d}\big{[}q_{a+e_{j}-e_{i},a}\,\tilde{\pi}_{n}(a+e_{j}-e_{i})-q_{a,a+e_{j}-e_{i}}\tilde{\pi}_{n}(a)\big{]}$
(3.10)
$\displaystyle+\nu(n)\sum_{i=1}^{d}\big{[}q_{a+e_{i},a}\frac{\sum_{i}\lambda_{i}}{\delta(n+1)}\tilde{\pi}_{n+1}(a+e_{i})-q_{a,a+e_{i}}\tilde{\pi}_{n}(a)\big{]}$
$\displaystyle+\nu(n)\sum_{i=1}^{d}\big{[}q_{a-e_{i},a}\frac{\delta\,n}{(\sum_{i}\lambda_{i})}\tilde{\pi}_{n-1}(a-e_{i})-q_{a,a-e_{i}}\tilde{\pi}_{n}(a)\big{]}$

where $\nu\sim$Poisson$(\sum_{i}\lambda_{i}/\delta)$ was used to express
$\frac{\nu(n+1)}{\nu(n)}$ and $\frac{\nu(n)}{\nu(n-1)}$ in the second and
third sums.

We now use
${\mathcal{B}}^{*}(a):={\mathcal{A}}^{*}\tilde{\Pi}(a)/\tilde{\Pi}(a)$ to
analyze how close $\tilde{\Pi}$ is to $\Pi$. By (3.10) we have

$\displaystyle{\mathcal{B}}^{*}(a)$ $\displaystyle=\sum_{i\neq
j=1}^{d}\big{[}q_{a+e_{j}-e_{i},a}\,\frac{\tilde{\pi}_{n}(a+e_{j}-e_{i})}{\tilde{\pi}_{n}(a)}-q_{a,a+e_{j}-e_{i}}\big{]}$
$\displaystyle+\sum_{i=1}^{d}\big{[}q_{a+e_{i},a}\frac{\sum_{i}\lambda_{i}}{\delta(n+1)}\frac{\tilde{\pi}_{n+1}(a+e_{i})}{\tilde{\pi}_{n}(a)}-q_{a,a+e_{i}}\big{]}+\sum_{i=1}^{d}\big{[}q_{a-e_{i},a}\frac{\delta\,n}{(\sum_{i}\lambda_{i})}\frac{\tilde{\pi}_{n-1}(a-e_{i})}{\tilde{\pi}_{n}(a)}-q_{a,a-e_{i}}\big{]}$

For $d=2$ each sum has two terms which we gather (based on $\pm$ sign) into

$\displaystyle{\mathcal{B}}^{*}(a)=\frac{1}{\tilde{\pi}_{n}(a)}\big{[}(L_{n-1}+L_{n}+L_{n+1})(a)-R_{n}(a)\big{]}$
with $\displaystyle
R_{n}(a)=[\lambda_{1}+\lambda_{2}+n\delta+(\kappa_{1}+\kappa_{2})a_{1}a_{2}]\tilde{\pi}_{n}(a)$
$\displaystyle
L_{n-1}(a)=\frac{n\delta\lambda_{1}}{\lambda_{1}+\lambda_{2}}\tilde{\pi}_{n-1}(a-e_{1})+\frac{n\delta\lambda_{2}}{\lambda_{1}+\lambda_{2}}\tilde{\pi}_{n-1}(a-e_{2})$
$\displaystyle
L_{n}(a)=\kappa_{2}(a_{1}+1)(a_{2}-1)\tilde{\pi}_{n}(a+e_{1}-e_{2})+\kappa_{1}(a_{1}-1)(a_{2}+1)\tilde{\pi}_{n}(a-e_{1}+e_{2})$
$\displaystyle
L_{n+1}(a)=\frac{\lambda_{1}+\lambda_{2}}{n+1}(a_{1}+1)\tilde{\pi}_{n+1}(a+e_{1})+\frac{\lambda_{1}+\lambda_{2}}{n+1}(a_{2}+1)\tilde{\pi}_{n+1}(a+e_{2})$

We see here the effect of the inflow and outflow reactions that give rise to
transitions from $n$ to $n\pm 1$, in different values of $n$ in $\tilde{\pi}$
and the Poisson distribution $\nu$. In the expression for $\tilde{\pi}$, the
weighted Dirichlet-multinomial depends only on $a\in E_{n}$, but the
normalizing factor $u$ depends on $n$, so we use Gauss hypergeometric
functions to simplify the ratios of $\tilde{\pi}$ appearing in
${\mathcal{B}}^{*}$.

In terms of a Dirichlet-multinomial$(n,\alpha)$ variable
$\xi=(\xi_{1},\dots,\xi_{d})$, the factor $u$ is, by (3.3)

$u(\alpha,\kappa,n)=\mathbb{E}\bigg{[}\bigg{(}\sum_{i=1}^{d}\kappa_{i}\xi_{i}\bigg{)}^{n}\bigg{]}$
(3.11)

which we can also write, using the product moments of $\xi$, as

$u(\alpha,\kappa,n)=\sum_{|a|=n}\bigg{[}\prod_{i=1}^{d}\kappa_{i}^{a_{i}}\bigg{]}{n\choose
a}\mathbb{E}\bigg{[}\prod_{i=1}^{d}\xi_{i}^{a_{i}}\bigg{]}=\sum_{|a|=n}\bigg{[}\prod_{i=1}^{d}\kappa_{i}^{a_{i}}\bigg{]}{n\choose
a}\frac{\Gamma(\sum_{i=1}^{d}\alpha_{i})}{\Gamma(\sum_{i=1}^{d}\alpha_{i}+n)}\prod_{i=1}^{d}\frac{\Gamma(\alpha_{i}+a_{i})}{\Gamma(\alpha_{i})}.$

For $d=2$ this can be further simplified using special functions. Using
$a_{1}=i,a_{2}=n-i$,

$u(\alpha,\kappa,n)=\sum_{i=0}^{n}{n\choose
i}\kappa_{1}^{i}\kappa_{2}^{n-i}\frac{\Gamma(\alpha_{1}+\alpha_{2})}{\Gamma(\alpha_{1}+\alpha_{2}+n)}\frac{\Gamma(\alpha_{1}+i)}{\Gamma(\alpha_{1})}\frac{\Gamma(\alpha_{2}+n-i)}{\Gamma(\alpha_{2})}$

and the partition function becomes

$\displaystyle u(\alpha,\kappa,n)$
$\displaystyle=\frac{\Gamma(\alpha_{1}+\alpha_{2})}{\Gamma(n+\alpha_{1}+\alpha_{2})}\frac{\Gamma(\alpha_{2}+n)}{\Gamma(\alpha_{2})}\kappa_{2}^{n}\sum_{i=0}^{n}{n\choose
i}(-1)^{i}\frac{(\alpha_{1})_{(i)}}{(1-\alpha_{2}-n)_{(i)}}\bigg{(}\frac{\kappa_{1}}{\kappa_{2}}\bigg{)}^{i}$
$\displaystyle=\frac{\Gamma(\alpha_{1}+\alpha_{2})}{\Gamma(n+\alpha_{1}+\alpha_{2})}\frac{\Gamma(\alpha_{2}+n)}{\Gamma(\alpha_{2})}\kappa_{2}^{n}{}_{2}F_{1}\big{(}-n,\alpha_{1};1-\alpha_{2}-n;\frac{\kappa_{1}}{\kappa_{2}}\big{)}$

where ${}_{2}F_{1}$ is the Gauss hypergeometric function, which when evaluated
on a nonpositive integer in the first coordinate reduces to the polynomial
${}_{2}F_{1}(-n,x;y;z)=\sum_{i=0}^{n}(-1)^{i}{n\choose
i}\frac{(x)_{(i)}}{(y)_{(i)}}{z^{i}}$.

The distribution $\tilde{\pi}_{n}$ on $i\in\\{0,1,\dots,n\\}$ becomes

$\displaystyle\tilde{\pi}_{n}(i)$
$\displaystyle=\frac{1}{u(\alpha,\kappa,n)}\kappa_{1}^{i}\kappa_{2}^{n-i}{n\choose
i}\frac{\Gamma(\alpha_{1}+\alpha_{2})}{\Gamma(\alpha_{1}+\alpha_{2}+n)}\frac{\Gamma(\alpha_{1}+i)\Gamma(\alpha_{2}+n-i)}{\Gamma(\alpha_{1})\Gamma(\alpha_{2})}$
$\displaystyle=\frac{1}{{}_{2}F_{1}\big{(}-n,\alpha_{1};1-\alpha_{2}-n;\frac{\kappa_{1}}{\kappa_{2}}\big{)}}\bigg{(}\frac{\kappa_{1}}{\kappa_{2}}\bigg{)}^{a}{n\choose
i}\frac{\Gamma(\alpha_{1}+i)\Gamma(\alpha_{2}+n-i)}{\Gamma(\alpha_{1})\Gamma(\alpha_{2}+n)}.$
(3.12)

Our representation by Gauss hypergeometric functions is valid when: (i) $y$ is
not a non-positive integer, and when (ii) $|z|<1$, or $|z|\leq 1$ when
$y>x-n$. This requires that: (i) $\alpha_{2}\notin\mathbb{N}-n$, and that (ii)
we label the species so that ${\kappa_{1}}/{\kappa_{2}}<1$, allowing for
$\kappa_{1}/\kappa_{2}\leq 1$ when $\delta<\kappa_{1}$. Since
$\alpha_{2}=(\delta\lambda_{2})/\kappa_{2}(\lambda_{1}+\lambda_{2})$ then
$\kappa_{1}/\kappa_{2}<1$ together with $\delta<\kappa_{1}$ will ensure that
$0<\alpha_{2}<1$ as well.

Using (3.1) in the expressions $R_{n},L_{n},L_{n-1},L_{n+1}$ for
${\mathcal{B}}^{*}$, we then get that (see Appendix A)

$\displaystyle{\mathcal{B}}^{*}(a)=$
$\displaystyle(\lambda_{1}+\lambda_{2})\Bigg{[}1-\frac{{}_{2}F_{1}(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{{}_{2}F_{1}(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}\Bigg{]}$
$\displaystyle-\kappa_{2}(n-1+\alpha_{2})\bigg{(}\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}+\frac{a_{2}\alpha_{2}}{a_{2}-1+\alpha_{2}}\bigg{)}\Bigg{[}\frac{{}_{2}F_{1}(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{{}_{2}F_{1}(1-n,\alpha_{1},2-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}-1\Bigg{]}$
$\displaystyle-\kappa_{2}\bigg{(}\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}+\frac{a_{2}\alpha_{2}}{a_{2}-1+\alpha_{2}}\bigg{)}\Big{[}a_{1}\big{(}1-\frac{\kappa_{1}}{\kappa_{2}})-\alpha_{2}\Big{]}-\kappa_{2}\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}\Big{[}\frac{\kappa_{1}}{\kappa_{2}}-1\Big{]}$
$\displaystyle-\frac{\lambda_{1}+\lambda_{2}}{n+\alpha_{2}}\frac{{}_{2}F_{1}(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{{}_{2}F_{1}(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}\Big{[}a_{1}\big{(}\frac{\kappa_{1}}{\kappa_{2}}-1\big{)}+\alpha_{2}\Big{]}$
(3.13)

We use numerical properties of Gauss hypergeometric functions ${}_{2}F_{1}$ to
determine the assumptions on parameters $\lambda_{i},\kappa_{i},$ and $\delta$
that ensure ${\mathcal{B}}^{*}(a)\approx 0$. We also investigate its
dependence on the overall species count in the system. Note that in the
symmetric case, $\kappa_{1}=\kappa_{2}$, simplifying
$R_{n},L_{n-1},L_{n},L_{n+1}$ yields ${\mathcal{B}}^{*}(a)=0$, $\forall a$,
and hence $\tilde{\Pi}$ is the true stationary distribution for $X$ (see proof
of Thm 4.2 in [20]).

###  3.2 Scaling and numerical results

We let $V$ denote the scaling parameter for the magnitude of the overall
species count in the system (referred to as volume in previous analyses of the
Togashi-Kaneko model). The classical mass-action scaling of the stochastic
model for the auto-catalytic reaction system (2.1) requires all bimolecular
reaction rate parameters to scale as $V^{-1}$, all inflow rate parameters to
scale as $V$, and all unimolecular outflow rate parameters to stay unscaled,

$\kappa_{i}=\frac{\kappa_{i}^{\prime}}{V},\;\;\lambda_{i}=\lambda_{i}^{\prime}V,\;\;\delta=\delta^{\prime}\;\;\quad\Rightarrow\quad\alpha_{i}=V\alpha_{i}^{\prime}$

To simplify the expression for
${\mathcal{B}}^{*}={\mathcal{A}}^{*}\tilde{\Pi}/\tilde{\Pi}$, and analyze its
properties for low versus high molecular count, we use the parametrization as
in [20], with $D,\kappa_{i}^{\prime}\sim O(1)$:

$\kappa_{i}=\frac{\kappa_{i}^{\prime}}{V},\;\;\lambda_{i}=DV,\;\;\delta=D,\;\;\quad\Rightarrow\quad\alpha_{i}=\frac{DV}{d\kappa_{i}^{\prime}}=V\alpha_{i}^{\prime}$
(3.14)

and focus on the role of asymmetric auto-catalytic rates $\kappa_{i}$ as well
as the effect of $V$. When $d=2$ we rewrite ${\mathcal{B}}^{*}$ from (3.1) in
terms of scaled parameters (see Appendix B) by factoring out the parameters
$DV$ and $\kappa_{1}^{\prime}/\kappa_{2}^{\prime}$ to get

$\begin{split}{\mathcal{B}}^{*}_{V}(a)=&dDV\Bigg{(}\Bigg{[}1-\frac{F(-n,\frac{DV}{d\kappa_{1}^{\prime}},1-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}{F(-1-n,\frac{DV}{d\kappa_{1}^{\prime}},-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}\Bigg{]}\\\
+&\frac{1}{n+\frac{DV}{d\kappa_{2}^{\prime}}}\frac{F(-n,\frac{DV}{d\kappa_{1}^{\prime}},1-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}{F(-1-n,\frac{DV}{d\kappa_{1}^{\prime}},-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}\Big{[}a_{1}\big{(}1-\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}}\big{)}-\frac{DV}{d\kappa_{2}^{\prime}}\Big{]}\Bigg{)}\\\
-&\kappa_{2}^{\prime}(n-1+\frac{DV}{d\kappa_{2}^{\prime}})\Bigg{[}\frac{F(-n,\frac{DV}{d\kappa_{1}^{\prime}},1-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}{F(1-n,\frac{DV}{d\kappa_{1}^{\prime}},2-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}-1\Bigg{]}\bigg{(}\frac{a_{1}\frac{DV}{d\kappa_{1}^{\prime}}}{a_{1}-1+\frac{DV}{d\kappa_{1}^{\prime}}}+\frac{a_{2}\frac{DV}{d\kappa_{2}^{\prime}}}{a_{2}-1+\frac{DV}{d\kappa_{2}^{\prime}}}\bigg{)}\\\
-&\kappa_{2}^{\prime}\Big{[}a_{1}\big{(}1-\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})-\frac{DV}{d}\Big{]}\bigg{(}\frac{a_{1}\frac{DV}{d\kappa_{1}^{\prime}}}{a_{1}-1+\frac{DV}{d\kappa_{1}^{\prime}}}+\frac{a_{2}\frac{DV}{d\kappa_{2}^{\prime}}}{a_{2}-1+\frac{DV}{d\kappa_{2}^{\prime}}}\bigg{)}+\kappa_{2}^{\prime}\frac{a_{1}\frac{D}{d\kappa_{1}^{\prime}}}{a_{1}-1+\frac{DV}{d\kappa_{1}^{\prime}}}\Big{[}1-\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}}\Big{]}.\end{split}$
(3.15)

First, note that in case
$\lambda_{i}^{\prime}=\delta^{\prime}=D\ll\kappa_{i}^{\prime}$ then
${\mathcal{B}}^{*}_{V}(a)=O(D)$ implies our proposed distribution
$\tilde{\Pi}$ is close to the stationary distribution, regardless of whether
the overall species count $V$ is low or high. Intuitively, when inflow and
outflow rates are small relative to the auto-catalytic rates, the system
spends most of its time with a constant overall species count, and hence the
distribution $\tilde{\pi}_{n}$ is stationary most of the time.

Figure 1: The ratios
$\frac{{}_{2}F_{1}(-n,\alpha_{1},1-\alpha_{2}-n,\frac{{\kappa_{1}}}{{\kappa_{2}}})}{{}_{2}F_{1}(1-n,\alpha_{1},2-\alpha_{2}-n,\frac{{\kappa_{1}}}{{\kappa_{2}}})}$
and
$\frac{{}_{2}F_{1}(-n,\alpha_{1},1-\alpha_{2}-n,\frac{{\kappa_{1}}}{{\kappa_{2}}})}{{}_{2}F_{1}(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{{\kappa_{1}}}{{\kappa_{2}}})}$
plotted in blue and orange, respectively, as a function of n. The rate
parameters are
$\lambda_{1}=\lambda_{2}=2,\,\delta=0.01,\,\kappa_{1}=1,\,\kappa_{2}=1.001$

Next, numerical explorations show the two ratios of the Gauss hypergeometric
functions appearing in ${\mathcal{B}}^{*}_{V}$ are both close to 1 for all
$(a_{1},a_{2})$ (Figure 1), so the two terms in the large square brackets are
small. Further, numerical evaluations of ${\mathcal{B}}^{*}_{V}(a)$ for
different values of $(a_{1},a_{2})$ are all close to zero (Figure 2). The very
small deviations from zero occur only when:
$(a_{1},a_{2})\in\\{(0,dV),(1,dV),(dV,0),(dV,1)\\}$. These are values with
$a_{1}+a_{2}\approx dV$ close to the mean of the overall species count
($\nu\sim$Poisson($\sum_{i}\lambda_{i}^{\prime}/\delta^{\prime}=dV)$).
Changing $\kappa^{\prime}_{1}/\kappa^{\prime}_{2}$ while keeping all other
parameters fixed, the deviations from zero disappear as
$\kappa^{\prime}_{1}/\kappa^{\prime}_{2}\to 1$ (see Figure 2) as expected.

Figure 2: $\mathcal{B}^{*}$ plotted with
$\lambda_{1}^{\prime}=\lambda_{2}^{\prime}=\delta^{\prime}=0.01,V=20$,
$\kappa_{1}^{\prime}=1$ and $\kappa_{2}^{\prime}=1.1,1.01,1.005,1.0001$
gradually decreasing $\kappa_{1}^{\prime}/\kappa_{2}^{\prime}\to 1$ from top
left to bottom right.

The effect of the scaling on $\tilde{\Pi}_{V}$ is reflected in both
$\nu_{V}\sim$Poisson($\sum_{i}\lambda_{i}^{\prime}/\delta^{\prime}=dV$) and in
the scaling of Gamma functions in $\tilde{\pi}_{n,V}$ (the $V$ scaling of
$\kappa_{i}^{\prime}$ in pre-factor cancels with same in $u_{V}$)

$\displaystyle\tilde{\pi}_{n,V}(a)=\frac{1}{u_{V}(\alpha^{\prime},\kappa^{\prime},n)}{n\choose
a}\Bigg{[}\prod_{i=1}^{d}(\kappa_{i}^{\prime})^{a_{i}}\Bigg{]}\frac{\Gamma\big{(}V\sum_{i=1}^{d}\alpha_{i}^{\prime}\big{)}}{\Gamma\big{(}n+V\sum_{i=1}^{d}\alpha_{i}^{\prime}\big{)}}\prod_{i=1}^{d}\frac{\Gamma\big{(}a_{i}+V\alpha_{i}^{\prime}\big{)}}{\Gamma\big{(}V\alpha_{i}^{\prime}\big{)}},$
$\displaystyle
u_{V}(\alpha^{\prime},\kappa^{\prime},n)=\frac{\Gamma\Big{(}\sum_{i=1}^{d}V\alpha_{i}^{\prime}\Big{)}}{\Gamma\Big{(}n+V\sum_{i=1}^{d}\alpha_{i}^{\prime}\Big{)}}\sum_{|a|=n}{n\choose
a}\prod_{i=1}^{d}(\kappa_{i}^{\prime})^{a_{i}}\frac{\Gamma(V\alpha_{i}^{\prime}+a_{i})}{\Gamma(V\alpha_{i}^{\prime})}$

For $d=2$ we can also write the scaled $\tilde{\pi}_{n,V}$ using $a=(i,n-i)$
in terms of Beta functions

$\tilde{\pi}_{n,V}(i)=\frac{1}{u_{V}(\alpha^{\prime},\kappa^{\prime},n)}(\kappa_{1}^{\prime})^{i}(\kappa_{2}^{\prime})^{(n-i)}{n\choose
i}\frac{B\big{(}i+\frac{DV}{d\kappa_{1}^{\prime}},n-i+\frac{DV}{d\kappa_{2}^{\prime}}\big{)}}{B\big{(}\frac{DV}{d\kappa_{1}^{\prime}},\frac{DV}{d\kappa_{2}^{\prime}}\big{)}}$
(3.16)

as a weighted Beta-binomial distribution. Just as in the symmetric case (see
Sec 4.3 [20]), our proposed distribution
$\tilde{\Pi}_{V}=\nu_{V}\tilde{\pi}_{n,V}$ will have a different shape:
unimodal or multimodal, depending on whether: $DV<d$ or $DV>d$.

When $DV<d$ the parameters
$\alpha^{\prime}_{1}=DV/(d\kappa_{1}^{\prime}),\alpha^{\prime}_{2}=DV/(d\kappa_{2}^{\prime})$
in the Beta-binomial are both less than 1, the weights have a small effect,
and the distribution is bimodal. Since
$\alpha^{\prime}_{2}/\alpha_{1}^{\prime}=\kappa_{1}^{\prime}/\kappa_{2}^{\prime}<1$,
the mode at the value $(a_{1},a_{2})=(0,n)$ is higher than at
$(a_{1},a_{2})=(n,0)$. Note that even very small asymmetries, e.g.
$\kappa_{2}^{\prime}=1.0001,\;\kappa^{\prime}_{1}=1$, show a sizeable
asymmetry in the size of the modes, and for a greater ratio, e.g.
$\kappa_{2}^{\prime}=1.1,\kappa^{\prime}_{1}=1$ almost the entire mass of the
distribution is concentrated at the higher mode (see Figure 3). Of course, on
each plane $(0,n)$ and $(n,0)$ there is a mode as well at $dV$ from the
distribution of total molecular counts $\nu_{V}$.

Figure 3: Plots of the analytic formula for $\tilde{\Pi}(a)$ when $DV=0.2<d=2$
show bimodal density. Parameters $D=0.01,V=20,\kappa_{1}^{\prime}=1$ are fixed
while $\kappa_{2}^{\prime}$ is increased in plots from left to right:
$\kappa_{2}^{\prime}=1.001,1.01,1.1$. Note the modes occur at
$(a_{1},a_{2})\in\\{(0,dV),(dV,0)\\}$ with $dV=40$.

Figure 4: Plots of $\tilde{\Pi}(a)$ when $DV=d=2$. Parameters are
$D=0.01,V=200,\kappa_{1}^{\prime}=1$, and $\kappa_{2}^{\prime}$ increases in
plots from left to right: $\kappa_{2}^{\prime}=1.001,1.01,1.1$.

In the intermediate case, when $DV=d=2$ we have $\alpha^{\prime}_{1}=1$ and
$\alpha_{2}^{\prime}=1/\kappa_{2}^{\prime}$. If $\kappa_{2}^{\prime}=1$ the
evaluation of $\tilde{\Pi}$ on hyperplanes $E_{n}$ would be Uniform($[0,n]$),
but for $\kappa_{2}^{\prime}>1$ it becomes skewed towards $\\{(0,n)\\}$, when
e.g. $\kappa_{2}=1.1$ almost all the mass is at $(0,dV=400)$, see Figure 4.

When $DV>2$ we have $\alpha^{\prime}_{1},\alpha^{\prime}_{2}>1$ and
$\tilde{\Pi}$ is unimodal concentrated around the fixed point
$(a_{1}^{*},a_{2}^{*})$ of the deterministic system for this model

$\begin{split}\frac{d}{dt}{a_{1}}(t)&=(\kappa^{\prime}_{1}-\kappa^{\prime}_{2})a_{1}a_{2}+\lambda^{\prime}_{1}-\delta^{\prime}a_{1}\\\
\frac{d}{dt}{a_{2}}(t)&=(\kappa^{\prime}_{2}-\kappa^{\prime}_{1})a_{1}a_{2}+\lambda^{\prime}_{2}-\delta^{\prime}a_{2}\end{split}$

With $\kappa_{1}^{\prime}=1$, increasing $\kappa_{2}^{\prime}>1$ the mode
moves such that the distribution is concentrated at a higher proportion of
$a_{2}$ species to $a_{1}$ species. When
$\kappa^{\prime}_{1}=\kappa^{\prime}_{2}$ this is a simple linear first order
ODE with fixed point
$(a_{1}^{*},a_{2}^{*})=(\lambda_{1}^{\prime}/\delta^{\prime},\lambda_{2}^{\prime}/\delta^{\prime})=(V,V).$
However, when $\kappa_{1}\neq\kappa_{2}$ this system of ODEs is nonlinear, and
we found numerically for $\kappa_{2}^{\prime}=1.001,1.01,1.1$ that the fixed
points are, respectively,
$(a_{1}^{*},a_{2}^{*})=(1801.96,2198.04),(763.93,3236.07),(97.5,3902.5)$
(Figure 5, left to right).

Figure 5: Plots of $\tilde{\Pi}(a)$ when $DV=20>d=2$. Parameters are
$D=0.01,V=2000,\kappa_{1}^{\prime}=1$, and $\kappa_{2}^{\prime}$ increases in
plots from left to right: $\kappa_{2}^{\prime}=1.001,1.01,1.1$.

Since the Markov chain $(X_{t})_{t\geq 0}$ is ergodic the long term fraction
of its occupation times converges to its stationary distribution. We simulated
the Markov chain, using the Gillespie algorithm, to further compare our
approximate stationary distribution $\tilde{\Pi}$ with the occupation time
when $DV<d$. In that case the modes at the boundaries of $E$ reflect the
Discreteness Induced Transition (DIT) nature of the system, where most of the
time is spent in the states $\\{(0,n),(n,0)\\}_{n\geq 0}$, and the quick
transitions between these states results in the lack of occupation time spent
in the interior of $E$. As the asymmetry
$\kappa_{2}^{\prime}/\kappa_{1}^{\prime}>1$ is increased, the DITs become more
and more rare, as the autocatalytic reaction
$A_{1}+A_{2}\xrightarrow{\kappa_{2}}2A_{2}$ overpowers the autocatalytic
reaction in the opposite direction. Figure 6 shows the empirical long-term
occupation time of the simulated Markov chain $X$ along with analytic plots of
$\tilde{\Pi}$.

On the other hand, DITs still exist but are less pronounced if we keep $DV<d$
but we increase the rates of inflows $\lambda_{i}^{\prime}$ and rate of
outflows $\delta^{\prime}$. In this case reactions
$\emptyset\xrightarrow{\lambda_{i}}A_{i}$, and
$A_{i}\xrightarrow{\delta}\emptyset$ compete with the effects of auto-
catalytic rates in either direction. The modes become less pronounced as the
process spends an equal amount of time in the interior of $E$ as in states
$\\{(0,n),(n,0)\\}_{n\geq 0}$. Figure 7 shows the effect of decreasing the
inflow and outflow rates while keeping the volume fixed, and shows the
empirical long term occupation time of the simulated Markov chain $X$ along
with analytic plots of $\tilde{\Pi}$. We note that despite the significant
effect of the inflow and outflow reactions, our distribution $\tilde{\Pi}$
matches the empirical estimates of $\Pi$ well on all of $E$. This agrees with
our numerical conclusions about $\mathcal{B}^{*}(a)\approx 0$ for all $a\in
E$.

Figure 6: Histograms for the fraction of occupation times in all states $a$ of
$X$ are shown in yellow with $\tilde{\Pi}(a)$ overlaid in blue. The parameters
are
$\lambda_{1}^{\prime}=\lambda_{2}^{\prime}=\delta^{\prime}=D=0.01,V=20,\kappa_{1}^{\prime}=1$,
varying $\kappa_{2}^{\prime}=1.001,1.01,1.1$ from left to right.

Figure 7: Histograms for the fraction of occupation times in all states $a$ of
$X$ are shown in yellow with $\tilde{\Pi}(a)$ overlaid in blue. The parameters
are
$\lambda_{1}^{\prime}=\lambda_{2}^{\prime}=D,\kappa_{1}^{\prime}=1,\kappa_{2}^{\prime}=1.01$,
and $V=20$ was fixed while $D$ was varied $D=2/100,5/100,7/100$ from left to
right.

##  4 Discussion

In this paper, we studied a family of auto-catalytic reaction systems
comprised of bimolecular interactions, plus inflows and outflows of molecular
species from the system. The stochastic Markov chain model for the system is
known to be ergodic with a unique stationary distribution, but the exact form
of this distribution is only known in the case of symmetric auto-catalytic
rates. To investigate the asymmetric case we proposed an approximate
stationary distribution $\tilde{\Pi}$ in an explicit form of a weighted
version of the distribution from the symmetric case. This distribution is
known to be the stationary distribution of a related model in population
genetics which is a conserved system with no inflows and outflows.

We derived an error function $\mathcal{B}^{*}$ which measures how close the
proposed distribution is to the true stationary one. Numerical evaluations of
this function, performed for a system with $d=2$ species, show that our
proposed distribution matches the true stationary distribution for almost all
values in the state space (and is very close to it in the exceptional values).
This is also confirmed by comparison of the empirical distribution of long-
term fraction of occupation times for the Markov chain with analytic values of
the proposed distribution.

The explicit form we provide for the proposed distribution $\tilde{\Pi}$
allows us to show that, under the usual mass-action scaling of the system
using the order of magnitude of the overall species count as a parameter, the
phenomenon of Discreteness Induced Transitions holds in the same parameter
regime as in the symmetric case. The asymmetry of the auto-catalytic rates
produces a significant asymmetry in the modes of the stationary distribution,
even under very small rate differences. The magnitude of inflow and outflow
rates relative to the auto-catalytic rates is shown to determine the extent to
which Discreteness Induced Transitions are felt in the system, with the speed
of transitions decreasing as the system becomes more open.

#### Acknowledgments

This work was initiated following an American Institute of Mathematics (AIM)
Workshop on Limits and Control of Stochastic Reaction Networks in August 2021,
and partial results were presented at the AIM Square group meeting in October
2022. This work constitutes the MSc thesis project of C.G. ans was supported
by L.P.’s NSERC Discovery Grant # RGPIN-2015-06573.

## References

  * [1] Bressloff, P. C. (2017). Stochastic switching in biology: from genotype to phenotype. Journal of Physics A: Mathematical and Theoretical, 50(13), 133001.
  * [2] Arcadia, C. E., Dombroski, A., Oakley, K., Chen, S. L., Tann, H., Rose, C., Eunsuk Kim, Sherief Reda, Brenda M. Rubenstein, & Rosenstein, J. K. (2021). Leveraging autocatalytic reactions for chemical domain image classification. Chemical Science, 12(15), 5464-5472.
  * [3] Hanopolskyi, A. I., Smaliak, V. A., Novichkov, A. I., & Semenov, S. N. (2021). Autocatalysis: Kinetics, mechanisms and design. ChemSystemsChem, 3(1), e2000026.
  * [4] Hordijk, W., & Steel, M. (2018). Autocatalytic networks at the basis of life’s origin and organization. Life, 8(4), 62.
  * [5] Bishop, L. M.,& Qian, H. (2010). Stochastic bistability and bifurcation in a mesoscopic signaling system with autocatalytic kinase. Biophysical journal, 98(1), 1-11.
  * [6] Vellela, M., & Qian, H. (2007). A quasistationary analysis of a stochastic chemical reaction: Keizer’s paradox. Bulletin of Mathematical Biology, 69, 1727-1746.
  * [7] Togashi, Y., & Kaneko, K. (2003). Alteration of chemical concentrations through discreteness-induced transitions in small autocatalytic systems. Journal of the Physical Society of Japan, 72(1), 62-68.
  * [8] Togashi, Y., & Kaneko, K. (2005). Discreteness-induced stochastic steady state in reaction diffusion systems: self-consistent analysis and stochastic simulations. Physica D: Nonlinear Phenomena, 205(1-4), 87-99.
  * [9] Togashi, Y., & Kaneko, K. (2007). Switching dynamics in reaction networks induced by molecular discreteness. Journal of Physics: Condensed Matter, 19(6), 065150.
  * [10] Biancalani, T., Rogers, T., & McKane, A. J. (2012). Noise-induced metastability in biochemical networks. Physical Review E, 86(1), 010106.
  * [11] Biancalani, T., Dyson, L., & Alan J McKane. (2014) Noise-induced bistable states and their mean switching time in foraging colonies. Physical review letters, 112(3):038101.
  * [12] Fan, W. T. L., Yang, Y. J., & Yuan, C. (2022). Constrained Langevin approximation for the Togashi-Kaneko model of autocatalytic reactions. arXiv preprint arXiv:2209.00426.
  * [13] Leite, S. & WIlliams, R. (2019) A constrained Langevin approximation for chemical reaction networks. The Annals of Applied Probability, 29(3):1541– 1608.
  * [14] K. Ball, T. Kurtz, , L. Popovic, & G. Rempala (2006). Multiscale Approximations for Reaction Networks. Annals of Applied Probability, 16 (4) pp 1925-1961.
  * [15] McSweeney, J. and Popovic, L. (2014). Stochastically-Induced Bistability in Chemical Reaction Systems, _Annals of Applied Probability_ , 24(3) pp 1226-1268.
  * [16] Hoessly, L. & Mazza, C. Stationary distributions and condensation in autocatalytic reaction networks. SIAM Journal on Applied Mathematics, 79(4):1173–1196, 2019.
  * [17] Nakagawa, M., & Togashi, Y. (2016). An analytical framework for studying small-number effects in catalytic reaction networks: A probability generating function approach to chemical master equations. Frontiers in Physiology, 7, 89.
  * [18] Yuichi Togashi and Kunihiko Kaneko. Transitions induced by the discreteness of molecules in a small autocatalytic system. Phys. Rev. Lett., 86:2459–2462, Mar 2001.
  * [19] Saito, N., & Kaneko, K. (2015). Theoretical analysis of discreteness-induced transition in autocatalytic reaction dynamics. Physical Review E, 91(2), 022707.
  * [20] E. Bibbona, J. Kim and C. Wiuf. Stationary distributions of systems with discreteness-induced transitions. J. R. Soc. Interface, 17, 2020.
  * [21] Etheridge, A.M., & Griffiths, R.C. (2009) A coalescent dual process in a Moran model with genic selection. Theoretical Population Biology, 75(4):320–330.
  * [22] Craciun, G., & Feinberg, M. (2006). Multiple equilibria in complex chemical reaction networks: extensions to entrapped species models. IEE Proceedings-Systems Biology, 153(4), 179-186.
  * [23] Banaji, M., & Pantea, C. (2018). The inheritance of nondegenerate multistationarity in chemical reaction networks. SIAM Journal on Applied Mathematics, 78(2), 1105-1130.
  * [24] Cappelletti, D., Feliu, E., & Wiuf, C. (2020). Addition of flow reactions preserving multistationarity and bistability. Mathematical biosciences, 320, 108295.

##  5 Appendices

###  5.1 Appendix A

We start from

$\displaystyle{\mathcal{B}}^{*}(a)=\frac{1}{\tilde{\pi}_{n}(a)}\big{[}(L_{n-1}+L_{n}+L_{n+1})(a)-R_{n}(a)\big{]}$
with $\displaystyle
R_{n}(a)=[\lambda_{1}+\lambda_{2}+n\delta+(\kappa_{1}+\kappa_{2})a_{1}a_{2}]\tilde{\pi}_{n}(a)$
$\displaystyle
L_{n-1}(a)=\frac{n\delta\lambda_{1}}{\lambda_{1}+\lambda_{2}}\tilde{\pi}_{n-1}(a-e_{1})+\frac{n\delta\lambda_{2}}{\lambda_{1}+\lambda_{2}}\tilde{\pi}_{n-1}(a-e_{2})$
$\displaystyle
L_{n}(a)=\kappa_{2}(a_{1}+1)(a_{2}-1)\tilde{\pi}_{n}(a+e_{1}-e_{2})+\kappa_{1}(a_{1}-1)(a_{2}+1)\tilde{\pi}_{n}(a-e_{1}+e_{2})$
$\displaystyle
L_{n+1}(a)=\frac{\lambda_{1}+\lambda_{2}}{n+1}(a_{1}+1)\tilde{\pi}_{n+1}(a+e_{1})+\frac{\lambda_{1}+\lambda_{2}}{n+1}(a_{2}+1)\tilde{\pi}_{n+1}(a+e_{2}).$

We use our expression (3.1) (dropping the subscripts in the Gauss
hypergeometric function ${}_{2}F_{1}$)

$\displaystyle\tilde{\pi}_{n}(i)$
$\displaystyle=\frac{1}{{}F\big{(}-n,\alpha_{1};1-\alpha_{2}-n;\frac{\kappa_{1}}{\kappa_{2}}\big{)}}\bigg{(}\frac{\kappa_{1}}{\kappa_{2}}\bigg{)}^{a}{n\choose
i}\frac{\Gamma(\alpha_{1}+i)\Gamma(\alpha_{2}+n-i)}{\Gamma(\alpha_{1})\Gamma(\alpha_{2}+n)}$

to get

$\displaystyle\frac{\tilde{\pi}_{n-1}(a-e_{1})}{\tilde{\pi}_{n}(a)}$
$\displaystyle=\frac{a_{1}}{n}\frac{n-1+\alpha_{2}}{a_{1}-1+\alpha_{1}}\frac{\kappa_{2}}{\kappa_{1}}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(1-n,\alpha_{1},2-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}$
$\displaystyle\frac{\tilde{\pi}_{n-1}(a-e_{2})}{\tilde{\pi}_{n}(a)}$
$\displaystyle=\frac{a_{2}}{n}\frac{n-1+\alpha_{2}}{a_{2}-1+\alpha_{2}}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(1-n,\alpha_{1},2-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}$
$\displaystyle\frac{\tilde{\pi}_{n}(a+e_{1}-e_{2})}{\tilde{\pi}_{n}(a)}$
$\displaystyle=\frac{a_{2}}{a_{1}+1}\frac{a_{1}+\alpha_{1}}{a_{2}-1+\alpha_{2}}\frac{\kappa_{1}}{\kappa_{2}}$
$\displaystyle\frac{\tilde{\pi}_{n}(a-e_{1}+e_{2})}{\tilde{\pi}_{n}(a+e_{1}-e_{2})}$
$\displaystyle=\frac{a_{1}}{a_{2}+1}\frac{a_{2}+\alpha_{2}}{a_{1}-1+\alpha_{1}}\frac{\kappa_{2}}{\kappa_{1}}$
$\displaystyle\frac{\tilde{\pi}_{n+1}(a+e_{1})}{\tilde{\pi}_{n}(a)}$
$\displaystyle=\frac{n+1}{a_{1}+1}\frac{a_{1}+\alpha_{1}}{n+\alpha_{2}}\frac{\kappa_{1}}{\kappa_{2}}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}$
$\displaystyle\frac{\tilde{\pi}_{n+1}(a)}{\tilde{\pi}_{n}(a)}$
$\displaystyle=\frac{n+1}{a_{2}+1}\frac{a_{2}+\alpha_{2}}{n+\alpha_{2}}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}$

Then $L_{n-1},L_{n},L_{n+1}$ divided by $\tilde{\pi}_{n}$ become

$\displaystyle\frac{L_{n-1}}{\tilde{\pi}_{n}(a)}=\frac{(n-1+\alpha_{2})\delta}{\lambda_{1}+\lambda_{2}}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(1-n,\alpha_{1},2-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}\Big{[}\frac{a_{1}\lambda_{1}}{a_{1}-1+\alpha_{1}}\frac{\kappa_{2}}{\kappa_{1}}+\frac{a_{2}\lambda_{2}}{a_{2}-1+\alpha_{2}}\Big{]}$
$\displaystyle\frac{L_{n}}{\tilde{\pi}_{n}(a)}=\kappa_{1}\frac{a_{2}-1}{a_{2}-1+\alpha_{2}}(a_{1}+\alpha_{1})a_{2}+\kappa_{2}\frac{a_{1}-1}{a_{1}-1+\alpha_{1}}(a_{2}+\alpha_{2})a_{1}$
$\displaystyle\frac{L_{n+1}}{\tilde{\pi}_{n}(a)}=\frac{\lambda_{1}+\lambda_{2}}{(n+\alpha_{2})}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}\Big{[}\frac{\kappa_{1}}{\kappa_{2}}(a_{1}+\alpha_{1})+(a_{2}+\alpha_{2})\Big{]}$

We use the fact that
$\alpha_{i}=\frac{\delta}{\lambda_{1}+\lambda_{2}}\frac{\lambda_{i}}{\kappa_{i}}$
(so $\sum_{i=1}^{2}\kappa_{i}\alpha_{i}=\delta$), $a_{1}+a_{2}=n$, and express
$\frac{L_{n-1}}{\tilde{\pi}_{n}}$ as

$\frac{L_{n-1}}{\tilde{\pi}_{n}(a)}=(n-1+\alpha_{2})\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(1-n,\alpha_{1},2-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}\kappa_{2}\Big{[}\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}+\frac{a_{2}\alpha_{2}}{a_{2}-1+\alpha_{2}}\Big{]}$
(5.1)

For $\frac{L_{n}}{\tilde{\pi}_{n}(a)}$ we arrange terms as follows

$\begin{split}\frac{L_{n}}{\tilde{\pi}_{n}(a)}=&\sum_{i=1}^{2}\sum_{j\neq
i}\kappa_{i}(a_{i}+\alpha_{i})\frac{a_{j}-1}{a_{j}-1+\alpha_{j}}a_{j}\\\
=&\sum_{i=1}^{2}\sum_{j\neq
i}\kappa_{i}(a_{i}+\alpha_{i})a_{j}-\sum_{i=1}^{2}\sum_{j\neq
i}\frac{\kappa_{i}(a_{i}+\alpha_{i})a_{j}\alpha_{j}}{a_{j}-1+\alpha_{j}}\\\
=&(\kappa_{1}+\kappa_{2})a_{1}a_{2}+\sum_{i=1}^{2}\sum_{j\neq
i}\kappa_{i}\alpha_{i}a_{j}-\sum_{i=1}^{2}\sum_{j\neq
i}\frac{\kappa_{i}(a_{i}+\alpha_{i})a_{j}\alpha_{j}}{a_{j}-1+\alpha_{j}}\\\
=&(\kappa_{1}+\kappa_{2})a_{1}a_{2}+\sum_{i=1}^{2}\sum_{j=1}^{2}\kappa_{i}\alpha_{i}a_{j}-\sum_{i=1}^{2}\kappa_{i}\alpha_{i}a_{i}-\sum_{i=1}^{2}\sum_{j=1}^{2}\frac{\kappa_{j}(a_{j}+\alpha_{j})a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}+\sum_{i=1}^{2}\frac{\kappa_{i}(a_{i}+\alpha_{i})a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\\\
=&(\kappa_{1}+\kappa_{2})a_{1}a_{2}+n\delta-\sum_{i=1}^{2}\kappa_{i}\alpha_{i}a_{i}-\sum_{i=1}^{2}\sum_{j=1}^{2}\frac{\kappa_{j}(a_{j}+\alpha_{j})a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}+\sum_{i=1}^{2}\kappa_{i}a_{i}\alpha_{i}+\sum_{i=1}^{2}\frac{\kappa_{i}a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\\\
=&(\kappa_{1}+\kappa_{2})a_{1}a_{2}+n\delta-\sum_{i=1}^{2}\sum_{j=1}^{2}\frac{\kappa_{j}(a_{j}+\alpha_{j})a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}+\sum_{i=1}^{2}\frac{\kappa_{i}a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\\\
=&(\kappa_{1}+\kappa_{2})a_{1}a_{2}+n\delta-\Big{(}\sum_{j=1}^{2}\kappa_{j}(a_{j}+\alpha_{j})\Big{)}\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}+\sum_{i=1}^{2}\frac{\kappa_{i}a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\\\
\end{split}$

so that the first two terms will cancel the last two terms in
$\frac{R_{n}}{\tilde{\pi}_{n}}$; for the remaining terms we have

$\begin{split}\quad&\sum_{i=1}^{2}\frac{\kappa_{i}a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}-\Big{(}\sum_{j=1}^{2}\kappa_{j}(a_{j}+\alpha_{j})\Big{)}\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\\\
&=\;\kappa_{1}\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}+\kappa_{2}\frac{a_{2}\alpha_{2}}{a_{2}-1+\alpha_{2}}\\\
&\;\;-\kappa_{1}(a_{1}+\alpha_{1})\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}-\kappa_{2}(a_{2}+\alpha_{2})\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}\\\
&\;\;-\kappa_{1}(a_{1}+\alpha_{1})\frac{a_{2}\alpha_{2}}{a_{2}-2+\alpha_{2}}-\kappa_{2}(a_{2}+\alpha_{2})\frac{a_{2}\alpha_{2}}{a_{2}-1+\alpha_{2}}\\\
&=\;\kappa_{2}\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}\Big{[}\frac{\alpha_{2}}{\alpha_{1}}-\frac{\alpha_{2}}{\alpha_{1}}(a_{1}+\alpha_{1})-(a_{2}+\alpha_{2})\Big{]}+\kappa_{2}\frac{a_{2}\alpha_{2}}{a_{2}-1+\alpha_{2}}\Big{[}1-\frac{\alpha_{2}}{\alpha_{1}}(a_{1}+\alpha_{1})-(a_{2}+\alpha_{2})\Big{]}\\\
&=\;\kappa_{2}\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}\Big{[}\frac{\alpha_{2}}{\alpha_{1}}-\frac{\alpha_{2}}{\alpha_{1}}a_{1}-a_{2}-2\alpha_{2}\Big{]}+\kappa_{2}\frac{a_{2}\alpha_{2}}{a_{2}-1+\alpha_{2}}\Big{[}1-\frac{\alpha_{2}}{\alpha_{1}}a_{1}-a_{2}-2\alpha_{2}\Big{]}\\\
&=\;\kappa_{2}\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{2}}\Big{[}\frac{\alpha_{2}}{\alpha_{1}}+a_{1}-\frac{\alpha_{2}}{\alpha_{1}}a_{1}-1-\alpha_{2}\Big{]}+\kappa_{2}\frac{a_{2}\alpha_{2}}{a_{2}-1+\alpha_{2}}\Big{[}a_{1}-\frac{\alpha_{2}}{\alpha_{1}}a_{1}-\alpha_{2}\Big{]}\\\
&\;\;-\kappa_{2}(n-1+\alpha_{2})\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\\\
&=\;\kappa_{2}\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\Big{[}a_{1}\big{(}1-\frac{\alpha_{2}}{\alpha_{1}}\big{)}-\alpha_{2}\Big{]}+\kappa_{2}\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}\Big{[}\frac{\alpha_{2}}{\alpha_{1}}-1\Big{]}-\kappa_{2}(n-1+\alpha_{2})\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\\\
&=\;\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\Big{[}a_{1}\big{(}\kappa_{2}-\kappa_{1})-\kappa_{2}\alpha_{2}\Big{]}+\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}\Big{[}\kappa_{1}-\kappa_{2}\Big{]}-\kappa_{2}(n-1+\alpha_{2})\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\\\
\end{split}$

so that

$\begin{split}\frac{L_{n}}{\tilde{\pi}_{n}(a)}=&(\kappa_{1}+\kappa_{2})a_{1}a_{2}+n\delta-\kappa_{2}(n-1+\alpha_{2})\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\\\
&+\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\Big{[}a_{1}\big{(}\kappa_{2}-\kappa_{1})-\alpha_{2}\Big{]}+\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}\Big{[}\kappa_{1}-\kappa_{2}\Big{]}\\\
\end{split}$ (5.2)

where the first two terms will cancel with terms in
$\frac{R_{n}}{\tilde{\pi}_{n}}$, and the third term is the same as in
$\frac{L_{n-1}}{\tilde{\pi}_{n}}$ from (5.1) apart from the ratio of
Hypergeometric functions.
Finally, we express $\frac{L_{n+1}}{\tilde{\pi}_{n}}$ as

$\begin{split}\frac{L_{n+1}}{\tilde{\pi}_{n}(a)}&=\frac{\lambda_{1}+\lambda_{2}}{n+\alpha_{2}}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}\Big{[}\frac{\kappa_{1}}{\kappa_{2}}(a_{1}+\alpha_{1})+(a_{2}+\alpha_{2})\Big{]}\\\
&=\frac{\lambda_{1}+\lambda_{2}}{n+\alpha_{2}}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}\Big{[}(n+\alpha_{2})+\frac{\kappa_{1}}{\kappa_{2}}(a_{1}+\alpha_{1})-a_{1}\Big{]}\\\
&=(\lambda_{1}+\lambda_{2})\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}+\frac{\lambda_{1}+\lambda_{2}}{n+\alpha_{2}}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}\Big{[}a_{1}\big{(}\frac{\kappa_{1}}{\kappa_{2}}-1\big{)}+\alpha_{2}\Big{]}\end{split}$

For ${\mathcal{B}}^{*}=\frac{1}{\tilde{\pi}_{n}}(L_{n-1}+L_{n}-L_{n-1}-R_{n})$
we now get the expression in (3.1)

$\begin{split}{\mathcal{B}}^{*}(a)=&(\lambda_{1}+\lambda_{2})\Bigg{[}1-\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}\Bigg{]}\\\
&-\kappa_{2}(n-1+\alpha_{2})\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\Bigg{[}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(1-n,\alpha_{1},2-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}-1\Bigg{]}\\\
&-\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}}{a_{i}-1+\alpha_{i}}\Big{[}a_{1}\big{(}\kappa_{2}-\kappa_{1})-\kappa_{2}\alpha_{2}\Big{]}-\frac{a_{1}\alpha_{1}}{a_{1}-1+\alpha_{1}}\Big{[}\kappa_{1}-\kappa_{2}\Big{]}\\\
&-\frac{\lambda_{1}+\lambda_{2}}{n+\alpha_{2}}\frac{F(-n,\alpha_{1},1-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}{F(-1-n,\alpha_{1},-\alpha_{2}-n,\frac{\kappa_{1}}{\kappa_{2}})}\Big{[}a_{1}\big{(}\frac{\kappa_{1}}{\kappa_{2}}-1\big{)}+\alpha_{2}\Big{]}.\end{split}$

###  5.2 Appendix B

We strat from (3.1) and use the (3.14) scaling by $V$

$\kappa_{i}=\frac{\kappa_{i}^{\prime}}{V},\;\;\lambda_{i}=DV,\;\;\delta=D,\;\;\quad\Rightarrow\quad\alpha_{i}=\frac{DV}{d\kappa_{i}^{\prime}}=V\alpha_{i}^{\prime}$

to get

$\begin{split}{\mathcal{B}}^{*}_{V}(a)=&dDV\Bigg{[}1-\frac{F(-n,V\alpha_{1}^{\prime},1-V\alpha_{2}^{\prime}-n,\frac{\frac{\kappa_{1}^{\prime}}{V}}{\frac{\kappa_{2}^{\prime}}{V}})}{F(-1-n,V\alpha_{1}^{\prime},-V\alpha_{2}^{\prime}-n,\frac{\frac{\kappa_{1}^{\prime}}{V}}{\frac{\kappa_{2}^{\prime}}{V}})}\Bigg{]}\\\
&-\frac{\kappa_{2}^{\prime}}{V}(n-1+V\alpha_{2}^{\prime})\sum_{i=1}^{2}\frac{a_{i}V\alpha_{i}^{\prime}}{a_{i}-1+V\alpha_{i}^{\prime}}\Bigg{[}\frac{F(-n,V\alpha_{1}^{\prime},1-V\alpha_{2}^{\prime}-n,\frac{\frac{\kappa_{1}^{\prime}}{V}}{\frac{\kappa_{2}^{\prime}}{V}})}{F(1-n,V\alpha_{1}^{\prime},2-V\alpha_{2}^{\prime}-n,\frac{\frac{\kappa_{1}^{\prime}}{V}}{\frac{\kappa_{2}^{\prime}}{V}})}-1\Bigg{]}\\\
&-\sum_{i=1}^{2}\frac{a_{i}V\alpha_{i}^{\prime}}{a_{i}-1+VV\alpha_{i}^{\prime\prime}}\Big{[}a_{1}\big{(}\frac{\kappa_{2}^{\prime}}{V}-\frac{\kappa_{1}^{\prime}}{V})-\frac{\kappa_{2}^{\prime}}{V}V\alpha_{2}^{\prime}\Big{]}-\frac{a_{1}V\alpha_{1}^{\prime}}{a_{1}-1+V\alpha_{1}^{\prime}}\Big{[}\frac{\kappa_{1}^{\prime}}{V}-\frac{\kappa_{2}^{\prime}}{V}\Big{]}\\\
&-\frac{\lambda_{1}+\lambda_{2}}{n+V\alpha_{2}^{\prime}}\frac{F(-n,V\alpha_{1}^{\prime},1-V\alpha_{2}^{\prime}-n,\frac{\frac{\kappa_{1}^{\prime}}{V}}{\frac{\kappa_{2}^{\prime}}{V}})}{F(-1-n,V\alpha_{1}^{\prime},-V\alpha_{2}^{\prime}-n,\frac{\frac{\kappa_{1}^{\prime}}{V}}{\frac{\kappa_{2}^{\prime}}{V}})}\Big{[}a_{1}\big{(}\frac{\frac{\kappa_{1}^{\prime}}{V}}{\frac{\kappa_{2}^{\prime}}{V}}-1\big{)}+V\alpha_{2}^{\prime}\Big{]}\\\
\end{split}$
$\begin{split}=&dDV\Bigg{[}1-\frac{F(-n,V\alpha_{1}^{\prime},1-V\alpha_{2}^{\prime}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}{F(-1-n,V\alpha_{1}^{\prime},-V\alpha_{2}^{\prime}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}\Bigg{]}\\\
&-\kappa_{2}(n-1+V\alpha_{2}^{\prime})\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}^{\prime}}{a_{i}-1+V\alpha_{i}^{\prime}}\Bigg{[}\frac{F(-n,V\alpha_{1}^{\prime},1-V\alpha_{2}^{\prime}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}{F(1-n,V\alpha_{1}^{\prime},2-V\alpha_{2}^{\prime}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}-1\Bigg{]}\\\
&-\sum_{i=1}^{2}\frac{a_{i}\alpha_{i}^{\prime}}{a_{i}-1+V\alpha_{i}^{\prime}}\Big{[}a_{1}\big{(}\kappa_{2}-\kappa_{1})-\kappa_{2}\alpha_{2}^{\prime}\Big{]}-\frac{a_{1}\alpha_{1}^{\prime}}{a_{1}-1+V\alpha_{1}^{\prime}}\Big{[}\kappa_{1}^{\prime}-\kappa_{2}^{\prime}\Big{]}\\\
&-\frac{dDV}{n+V\alpha_{2}^{\prime}}\frac{F(-n,V\alpha_{1}^{\prime},1-V\alpha_{2}^{\prime}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}{F(-1-n,V\alpha_{1}^{\prime},-V\alpha_{2}^{\prime}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}\Big{[}a_{1}\big{(}\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}}-1\big{)}+V\alpha_{2}^{\prime}\Big{]}\\\
\end{split}$

Substituting $\alpha_{i}^{\prime}=\frac{D}{d\kappa_{i}^{\prime}}$ and
rearranging terms we get the expression in (3.15)

$\begin{split}{\mathcal{B}}^{*}_{V}(a)=&dDV\Bigg{(}\Bigg{[}1-\frac{F(-n,\frac{DV}{d\kappa_{1}^{\prime}},1-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}{F(-1-n,\frac{DV}{d\kappa_{1}^{\prime}},-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}\Bigg{]}\\\
+&\frac{1}{n+\frac{DV}{d\kappa_{2}^{\prime}}}\frac{F(-n,\frac{DV}{d\kappa_{1}^{\prime}},1-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}{F(-1-n,\frac{DV}{d\kappa_{1}^{\prime}},-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}\Big{[}a_{1}\big{(}1-\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}}\big{)}-\frac{DV}{d\kappa_{2}^{\prime}}\Big{]}\Bigg{)}\\\
-&\kappa_{2}^{\prime}(n-1+\frac{DV}{d\kappa_{2}^{\prime}})\Bigg{[}\frac{F(-n,\frac{DV}{d\kappa_{1}^{\prime}},1-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}{F(1-n,\frac{DV}{d\kappa_{1}^{\prime}},2-\frac{DV}{d\kappa_{2}^{\prime}}-n,\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})}-1\Bigg{]}\bigg{(}\frac{a_{1}\frac{DV}{d\kappa_{1}^{\prime}}}{a_{1}-1+\frac{DV}{d\kappa_{1}^{\prime}}}+\frac{a_{2}\frac{DV}{d\kappa_{2}^{\prime}}}{a_{2}-1+\frac{DV}{d\kappa_{2}^{\prime}}}\bigg{)}\\\
-&\kappa_{2}^{\prime}\Big{[}a_{1}\big{(}1-\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}})-\frac{DV}{d}\Big{]}\bigg{(}\frac{a_{1}\frac{DV}{d\kappa_{1}^{\prime}}}{a_{1}-1+\frac{DV}{d\kappa_{1}^{\prime}}}+\frac{a_{2}\frac{DV}{d\kappa_{2}^{\prime}}}{a_{2}-1+\frac{DV}{d\kappa_{2}^{\prime}}}\bigg{)}+\kappa_{2}^{\prime}\frac{a_{1}\frac{D}{d\kappa_{1}^{\prime}}}{a_{1}-1+\frac{DV}{d\kappa_{1}^{\prime}}}\Big{[}1-\frac{\kappa_{1}^{\prime}}{\kappa_{2}^{\prime}}\Big{]}.\end{split}$