# Synchronization and random attractors for reaction jump processes

Maximilian<EMAIL_ADDRESS>Freie Universität Berlin,
Institut für Mathematik und Informatik, 14195 Berlin, Germany  Guillermo
Olicón-Méndez  Freie Universität Berlin, Institut für Mathematik und
Informatik, 14195 Berlin, Germany  Nathalie Unger  Zuse Institute Berlin,
14195 Berlin, Germany  Technische Universität Berlin, 10623 Berlin, Germany
Stefanie<EMAIL_ADDRESS>Zuse Institute Berlin, 14195 Berlin,
Germany

###### Abstract

This work explores a synchronization-like phenomenon induced by common noise
for continuous-time Markov jump processes given by chemical reaction networks.
A corresponding random dynamical system is formulated in a two-step procedure,
at first for the states of the embedded discrete-time Markov chain and then
for the augmented Markov chain including also random jump times. We uncover a
time-shifted synchronization in the sense that – after some initial waiting
time – one trajectory exactly replicates another one with a certain time
delay. Whether or not such a synchronization behaviour occurs depends on the
combination of the initial states. We prove this partial time-shifted
synchronization for the special setting of a birth-death process by analyzing
the corresponding two-point motion of the embedded Markov chain and determine
the structure of the associated random attractor. In this context, we also
provide general results on existence and form of random attractors for
discrete-time, discrete-space random dynamical systems.

Keywords: chemical reaction networks, random attractors, random periodic
orbits, reaction jump processes, synchronization

MSC2020: 37H99, 60J27, 60J10, 92C40

##  1 Introduction

Stochastic models of biochemical reaction dynamics are mostly based on the
theory of Markov processes [1, 26]. A central role play reaction jump
processes which model a well-mixed reaction system as a continuous-time Markov
process on a discrete state space. The state is given by the number of
particles of each involved chemical species, and chemical reactions are
modeled as stochastic events which induce jumps in the system’s state that
occur after exponentially distributed sojourn times. The temporal evolution of
the system’s probability distribution is in this case characterized by the
well-known chemical master equation [14]. Besides such reaction jump
processes, there exist also modeling approaches using discrete-time Markov
chains [19, 16], stochastic differential equations (SDEs) or ordinary
differential equations (ODEs) [15, 22, 21], which approximate the dynamics on
a macroscopic level in case of large population sizes, as well as hybrid model
recombinations for multiscale reaction systems [17, 37, 33, 27]. All these
approaches for describing and analyzing stochastic phenomena within
biochemical or other types of applied contexts have extensively been studied
in the literature [34, 2, 32].

###  1.1 Background and related work

The counterpart to stochastic processes within dynamical system theory is
given by random dynamical systems (RDS). Here, the origin of uncertainties is
considered somewhat differently. In simple terms, the system evolves according
to deterministic maps which are chosen randomly from a stochastic law.
Formally speaking, an RDS $(\theta,\varphi)$ on a metric state space
$\mathbb{X}$ (endowed with its Borel $\sigma$-algebra
$\mathcal{B}(\mathbb{X})$) and discrete time set $\mathbb{T}=\mathbb{N}_{0}$
or $\mathbb{Z}$ consists of

  * •

a noise model, given by a metric dynamical system
$(\Omega,\mathcal{F},\mathbb{P},(\theta_{n})_{n\in\mathbb{T}})$. By this we
mean that $(\Omega,\mathcal{F},\mathbb{P})$ is a probability space and
$\theta:=(\theta^{n})_{n\in\mathbb{T}}$ is a family of measurable maps
$\theta^{n}:\Omega\rightarrow\Omega$ for which
$\theta^{n+m}=\theta^{n}\circ\theta^{m}$ for all $n,m\in\mathbb{T}$, and which
is invariant with respect to $\mathbb{P}$ (or, $\mathbb{P}$ is
$\theta$-invariant). This last statement means that
$\theta^{n}_{*}\mathbb{P}(\cdot):=\mathbb{P}((\theta^{n})^{-1}(\cdot))=\theta(\cdot)$,

  * •

a cocycle map
$\varphi:\mathbb{N}_{0}\times\Omega\times\mathbb{X}\rightarrow\mathbb{X}$,
with $(n,\omega,x)\mapsto\varphi^{n}_{\omega}(x)$, which is measurable and
satisfies the cocycle property over $\theta$, that is for all
$x\in\mathbb{X}$, $n,m\in\mathbb{N}_{0}$, and $\omega\in\Omega$

$\varphi^{0}_{\omega}(x)=x,\qquad\varphi^{n+m}_{\omega}(x)=\varphi^{n}_{\theta^{m}\omega}\circ\varphi^{m}_{\omega}(x).$
(1)

Notice that while the dynamics on the noise space $\Omega$ might be defined
for both positive and negative times, this does not need to be the case for
the cocycle map $\varphi$ which in our context will only be defined for times
on $\mathbb{N}_{0}$. For a comprehensive theoretical background of RDS we
refer to [3].

For a discrete-time system on a finite state space the maps are given by
deterministic transitions matrices (containing only entries zero and one), and
the expectation of the matrix-valued random variable of transitions maps
agrees with the stochastic transition matrix of the corresponding Markov
chain. The relation between such finite-state RDS and the related Markov
chains has been studied by F. Ye et al. [36, 35]. Among other things, it has
been found that a given finite-state RDS induces a unique Markov chain, while
one Markov chain might be compatible with several RDS [36], as already
discussed in a general context by Kifer [18]. In this sense, the RDS
formulation may be seen as a more refined model of stochastic dynamics than
the Markov chain: the former gives a precise description of the two-point
motion, comparing trajectories with different initial conditions but driven by
the same noise allowing for the analysis of _random attractors_ [7], whereas
the latter characterizes the statistics of the one-point motion by means of
the transition probabilities.

While RDS representations of Markov chains (discrete in space and time) or
SDEs (continuous in space and time) (see e.g. [3]) have been studied in the
literature, an analogous investigation for continuous-time Markov processes on
discrete state spaces is still missing. In the present work, we do a first
step in this direction by formulating random dynamical systems corresponding
to reaction jump processes as special types of continuous-time Markov
processes. Our goal is to study questions of synchronization: Given the same
noise realization, will trajectories starting at different initial states
approach each other in the course of time? Once they coincide at a certain
time point, do they stay together forever? Numerical experiments have shown
that two realizations of the reaction jump process with distinct starting
points (but the same driving noise) may actually resemble each other after
some time period in the sense that one of the trajectories appears to be a
time-delayed replicate of the other. That is, after some random initial
“finding time”, the two process realizations start to wander through the same
sequence of states, with identical sojourn times in each of these states, but
with a certain time lag with respect to each other. Whether or not this type
of trajectory replication happens seems to depend in general on the
combination of chosen initial states. By means of the RDS presentation of the
dynamics, we provide an analytical explanation for this intriguing phenomenon
of time-shifted synchronization and its dependency on the initial conditions.

###  1.2 Main results

For our analysis, we use the fact that a (continuous-time) Markov jump process
$(X(t))_{t\geq 0}$ has a discrete-time representation given by the augmented
Markov chain [31] which assigns to each discrete index $n$ the random time
$T_{n}$ where the $n$th jump of the process occurs, as well as the state
$X_{n}=X(T_{n})$ entered by the process at this jump time. The random sequence
$(X_{n})_{n\in\mathbb{N}_{0}}$ of states, called embedded Markov chain, is a
standard (discrete-time) Markov chain on a countable state space. In
particular, we use an explicit recursive formula for this Markov chain which
immediately yields the cocycle of an RDS. We show that the time-shifted
synchronization of a time-homogeneous reaction jump process is equivalent to
the “normal” synchronization of the embedded Markov chain, for an appropriate
subset of initial conditions. Given that the jump rate constants are time-
independent, also the sojourn times within the states will agree once that the
states do agree.

In more detail, we focus on two main examples, providing several general
insights on random attractors for discrete state spaces on the way: a simple
birth-death process and the Schlögl model with their monostable and bistable
structures, respectively, detecting similarities and differences in the
described synchronization behaviour. We obtain the following main results and
insights:

  * •

For the embedded Markov chain of the birth-death process we prove _partial
synchronization_ (and, by that, partial time-shifted synchronization for the
reaction jump process) in the sense that common-noise trajectories with
starting states of the same parity join each other in finite time while
initial states with different parity lead to oscillations around each other
(Proposition 2 and Corollary 3).

  * •

For general RDS corresponding with Markov chains, we relate different forms of
random attractors (Theorem 8) and give conditions for the existence of a
(weak) random attractor (Theorem 10). We verify these conditions for our
examples (Proposition 11), prominently using the existence of a unique
stationary distribution.

  * •

We provide a main analytical result for the birth-death case, characterizing
the random attractor as a pullback and forward attractor consisting of two
random points with distance $1$ that form a random periodic orbit (Theorem
14), by that also finding the structure of the corresponding sample measures,
also called _statistical equilibria_ (Proposition 15).

  * •

We illustrate numerical insights that the weak attractor for the Schlögl model
has the same structure as for simple birth-death, apart from the fact that the
distance of the two random points is not $1$, mirroring the bistability of the
model.

Except for the general insights on random attractors on countable state
spaces, most of our analytical results are, so far, restricted to the special
case of a birth-death process since the absorbing property of the (thick)
diagonal for the two-point motion can be used for this case. However, the
general structure of the proof may well be extended to more general chemical
reaction networks, also with multiple reactants and corresponding random
periodic orbits.

Note that the works [36, 35, 16] mentioned earlier also deal with
synchronization of RDSs for Markov chains, and in [16] even partial
synchronization is considered. However, the latter approach is focused on
linear cocycles for random networks, using the theory of Lyapunov exponents.
Our proofs deploy an analysis of the two-point motion and its consequences for
the random attractor, and do not require a linear interpretation, being
confronted with an infinite state space. Note that Newman’s work on
synchronization for RDS [28, 29] achieves general equivalent conditions for
synchronization to occur, which can typically be verified via the maximal
Lyapunov exponent when the state space is a smooth manifold. For the class of
examples considered in this work, the equivalence of these conditions,
adjusted to the problem of partial synchronization, will automatically appear
in a straight-forward manner.

###  1.3 Structure of the paper

The remainder of the paper is structured as follows. In Sec. 2, we introduce
reaction jump processes and the corresponding augmented and embedded Markov
chains, interpreting the latter as random dynamical systems, and relate their
different forms of (partial) synchronization. Sec. 3 is dedicated to proving
partial synchronization of the birth-death chain, by using arguments on the
two-point motion. In Sec. 4, we discuss general properties of weak, pullback
and forward attractors for the discrete setting (Sec. 4.1), show a general
result on the existence of weak attractors including the setting of reaction
problems (Sec. 4.2), characterize the structure of this attractor for the
birth-death case as a random periodic orbit (Sec. 4.3) and give illustrations
of the two-point motion and the stationary, statistical behavior also for the
more complicated Schlögl model (Sec. 4.4). Finally, we provide a conclusion
with outlook in Sec. 5.

##  2 Reaction jump processes and random dynamical systems

In this section, we introduce the reaction network under consideration and
formulate the corresponding stochastic dynamics. At first (in Sec. 2.1), the
pathwise formulation of the reaction jump process is given, including two
exemplary reaction networks which will be extensively studied in this work.
The related random dynamical systems will be formulated in Sec. 2.2 and Sec.
2.3.

###  2.1 The reaction jump process

We consider the standard setting of well-mixed stochastic chemical reaction
dynamics [34]: There is a system of particles with $L\in\mathbb{N}$ different
types/species $\mathcal{S}_{1},\dotsc,\mathcal{S}_{L}$. The particles interact
by $K\in\mathbb{N}$ chemical reactions
$\mathcal{R}_{1},\dotsc,\mathcal{R}_{K}$ given by

$\mathcal{R}_{k}:\quad\sum_{l=1}^{L}s_{lk}\mathcal{S}_{l}\to\sum_{l=1}^{L}s^{\prime}_{lk}\mathcal{S}_{l},$

where the stoichiometric coefficients $s_{kl},s^{\prime}_{lk}$ are non-
negative integers. The state of the system is given by a vector
$x=(x_{l})_{l=1,...,L}\in\mathbb{N}_{0}^{L}$ with $x_{l}$ counting the number
of particles of species $l$. Each reaction induces a jump in the state of the
form $x\mapsto x+\nu_{k}$, where
$\nu_{k}=(\nu_{1k},...,\nu_{Lk})\in\mathbb{Z}^{L}$ is the state-change vector
given by $\nu_{lk}:=s^{\prime}_{kl}-s_{lk}$. Given a state $x$, the reaction
$\mathcal{R}_{k}$ takes place at rate $\alpha_{k}(x)$, where
$\alpha_{k}:\mathbb{N}_{0}^{L}\to[0,\infty)$ is the corresponding propensity
function.

The resulting reaction jump process (RJP) on $\mathbb{X}=\mathbb{N}_{0}^{L}$
has the path-wise representation

$X(t)=X(0)+\sum_{k=1}^{K}\mathcal{U}_{k}\left(\int_{0}^{t}\alpha_{k}(X(s))ds\right)\nu_{k},$

where $\mathcal{U}_{k}$ are independent unit-rate Poisson processes. This
process (and equivalently a more general Markov jump process) is fully
characterized by the random jump times $T_{n}$, $n=1,2,...$, at which the
jumps (here reactions) take place and the states $X_{n}:=X(T_{n})$ that are
entered at the jump times, namely by

$X(t)=X_{n}\quad\mbox{for}\;T_{n}\leq t<T_{n+1},$

with $T_{0}=0$ and $X_{0}=X(0)$. That is, we can consider a division of the
Markov jump process into the process of jump times
$(T_{n})_{n\in\mathbb{N}_{0}}$ with values in $[0,\infty)$ and the process of
the states $(X_{n})_{n\in\mathbb{N}_{0}}$ in $\mathbb{X}$ which is called the
embedded Markov chain. The discrete-time process
$(X_{n},T_{n})_{n\in\mathbb{N}_{0}}$ is called the augmented Markov chain
[31].

The jump times and the states of the RJP are recursively given by

$\displaystyle T_{n+1}$ $\displaystyle=T_{n}+\tau(X_{n})$ (2) $\displaystyle
X_{n+1}$ $\displaystyle=X_{n}+\nu_{\kappa(X_{n})},$ (3)

with $T_{0}=0$, $X_{0}=X(0)$, where $\tau(x)$ is an exponentially distributed
random variable with mean $1/x$ and $\kappa(x)\in\\{1,...,K\\}$ is a random
variable with point probabilities $\alpha_{k}(x)/\sum_{l=1}^{K}\alpha_{l}(x)$
for $k=1,...,K$. It is well known (see [12, 13]) that $\tau$ and $\kappa$ can
be realised by taking independent, uniformly distributed random numbers
$r,q\sim U(0,1)$ and setting

$\tau(x,r)=\frac{1}{\sum_{k=1}^{K}\alpha_{k}(x)}\log\left(\frac{1}{r}\right),$
(4)

and $\kappa(x,q)$ is the smallest integer satisfying

$\sum_{k=1}^{\kappa(x,q)}\alpha_{k}(x)>q\sum_{k=1}^{K}\alpha_{k}(x).$ (5)

Note that on a pathwise level $\kappa(x,q)$ depends on the order of reaction
indices.

######  Example 1 (Birth-death process).

As a basic example which will be analyzed in detail in Sec. 3.2 and Sec. 4.3
we consider the standard birth-death process of a single species $\mathcal{S}$
given by $K=2$ reactions

$\mathcal{R}_{1}:\emptyset\stackrel{{\scriptstyle\gamma_{1}}}{{\longrightarrow}}\mathcal{S},\quad\mathcal{R}_{2}:\mathcal{S}\stackrel{{\scriptstyle\gamma_{2}}}{{\longrightarrow}}\emptyset.$

Here, $\gamma_{1},\gamma_{2}>0$ are rate constants and the corresponding
propensity functions are given by the law of mass action as

$\alpha_{1}(x)=\gamma_{1},\quad\alpha_{2}(x)=\gamma_{2}x.$

The state space of the resulting jump process is given by
$\mathbb{X}=\mathbb{N}_{0}$. Consequently, also the state-change vectors
$\nu_{k}$ are actually scalar and given by $\nu_{1}=1$ and $\nu_{2}=-1$. From
(5) we can deduce that

$\kappa(x,q)=\begin{cases}1&\textup{if
}q<\frac{\gamma_{1}}{\gamma_{1}+\gamma_{2}x},\\\
2&\textup{otherwise,}\end{cases}$ (6)

which means that reaction $\mathcal{R}_{1}$ takes place with probability
$\frac{\gamma_{1}}{\gamma_{1}+\gamma_{2}x}$ given that the system is in state
$x$, while $\mathcal{R}_{2}$ takes place with probability
$1-\frac{\gamma_{1}}{\gamma_{1}+\gamma_{2}x}$.

Under an appropriate scaling of the propensity functions $\alpha_{1}$ and
$\alpha_{2}$, one may derive the corresponding _reaction rate equation_
governing the dynamics of the concentration $C(t)=\frac{X^{V}(t)}{V}$ for the
large volume limit $V\to\infty$ (cf. [34]). This is an ordinary differential
equation (ODE), given by

$\frac{\mathrm{d}C(t)}{\mathrm{d}t}=-\gamma_{2}C(t)+\gamma_{1},$ (7)

with globally attracting equilibrium at $C=\frac{\gamma_{1}}{\gamma_{2}}$ (see
Figure 1).

######  Example 2 (Schlögl model).

The other main example of this work is the Schlögl model, a chemical reaction
network exhibiting bistability, cf. e.g. [9, 30, 25] and see Sec. 4.4 for a
detailed discussion. Again we only have one species $\mathcal{S}$ with the
following reactions

$\begin{array}[]{lcl}\mathcal{R}_{1}:\emptyset\stackrel{{\scriptstyle\gamma_{1}}}{{\longrightarrow}}\mathcal{S},&&\mathcal{R}_{2}:\mathcal{S}\stackrel{{\scriptstyle\gamma_{2}}}{{\longrightarrow}}\emptyset,\\\
\mathcal{R}_{3}:2\mathcal{S}\stackrel{{\scriptstyle\gamma_{3}}}{{\longrightarrow}}3\mathcal{S},&&\mathcal{R}_{4}:3\mathcal{S}\stackrel{{\scriptstyle\gamma_{4}}}{{\longrightarrow}}2\mathcal{S},\end{array}$

for rate constants $\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}>0$ and
$x\in\mathbb{N}_{0}$. The corresponding standard mass action propensities are
given by

$\begin{array}[]{lcl}\alpha_{1}(x)=\gamma_{1},&&\alpha_{2}(x)=\gamma_{2}x,,\\\
\alpha_{3}(x)=\gamma_{3}x(x-1),&&\alpha_{4}(x)=\gamma_{4}x(x-1)(x-2).\end{array}$

The Schlögl model has the state-change vectors $\nu_{1},\nu_{3}=1$ and
$\nu_{2},\nu_{4}=-1$.

The corresponding reaction rate ODE is given by

$\frac{\mathrm{d}C(t)}{\mathrm{d}t}=-\gamma_{4}C(t)^{3}+\gamma_{3}C(t)^{2}-\gamma_{2}C(t)+\gamma_{1},$
(8)

which – depending on the value of the reaction rates – may exhibit an unstable
equilibrium and two stable equilibria (see Figure 1).

The two examples demonstrate two fundamentally different patterns in terms of
the large-volume behavior of the process: whereas in Example 1 the reaction
rate equation has one globally attracting equilibrium such that all
trajectories synchronize to the same concentration, in Example 2 there are two
(locally) attracting equilibria separated by an unstable equilibrium such that
different ODE-trajectories may or may not synchronize depending on their
initial conditions, see Figure 1.

This observation prompts interest in the synchronizing behavior of the
underlying reaction jump process: Given the same random numbers but different
initial states, will the jump times and states given by (2)-(3) approach each
other in the course of time? Can we observe different synchronization behavior
in Examples 1 and 2? To give a systematic answer to these questions, we will
analyze the reaction jump process in terms of the corresponding RDS which
gives a natural approach to comparing trajectories with different initial
conditions but driven by the same noise realizations.

(a)

(b)

Figure 1: Solution of the reaction rate equation for (a) the birth-death
process, see Eq. (7), and (b) the Schlögl model, see Eq. (8), each for three
different initial states $C(0)$. (a) Rate constants
$\gamma_{1}=10,\gamma_{2}=1$. (b) Rate constants
$\gamma_{1}=6,\gamma_{2}=3.5,\gamma_{3}=0.4,\gamma_{4}=0.0105$. Unstable
equilibrium given by $c\approx 9.6201$.

###  2.2 RDS formulation of the embedded Markov chain

At first, we formulate the setting of a RDS for the embedded Markov chain
$(X_{n})_{n\geq 0}$, given by $X_{n}:=X(T_{n})$, as a discrete-time stochastic
process.

The noise space $\mathcal{Q}_{+}$ of the RDS is chosen as

$\mathcal{Q}_{+}=\\{q=(q_{n})_{n\in\mathbb{N}_{0}}:q_{n}\in[0,1]\\}.$

We endow $\mathcal{Q}_{+}$ with the Borel $\sigma$-algebra
$\sigma(\mathcal{Q}_{+})$ generated by its cylinder sets, and with the
infinite product probability measure $\mathbb{P}=\lambda^{\mathbb{N}_{0}}$,
where $\lambda$ denotes the Lebesgue measure on $[0,1]$. On this probability
space $(\mathcal{Q}_{+},\sigma(\mathcal{Q}_{+}),\mathbb{P})$ we define the
shift map $\theta:\mathcal{Q}_{+}\to\mathcal{Q}_{+}$ and its iterates by

$\theta(q_{0},q_{1},...)=(q_{1},q_{2},...),\qquad\theta^{n}:=\underbrace{\theta\circ\cdots\circ\theta}_{n\textup{
times }}$ (9)

Since $\theta$ is invariant with respect to $\mathbb{P}$, the tuple
$(\mathcal{Q}_{+},\sigma(\mathcal{Q}_{+}),\mathbb{P},(\theta)_{n\in\mathbb{N}_{0}})$
constitutes our underlying noise model. Throughout this work we will use
interchangeably the short-hand notation

$\mathbb{P}(S(q))=\mathbb{P}(\\{q\in\mathcal{Q}_{+}:S(q)\textup{ holds}\\}),$

where $S(q)$ is a $q$-dependent statement. For any $q\in\mathcal{Q}_{+}$,
consider the transition map $f_{q}:\mathbb{X}\to\mathbb{X}$ defined by

$f_{q}(x):=x+\nu_{\kappa(x,q_{0})},$ (10)

where $f_{\cdot}$ takes a whole sequence $q=(q_{n})_{n\in\mathbb{N}_{0}}$ as
an input but only evaluates the first entry of $q$, namely $q_{0}$, in
$\kappa(x,\cdot)$ defined in (5). Therefore $f_{\theta^{n}q}(X_{n})$ coincides
with the right-hand side of the recursion (3). Given a fixed order of the
reactions, the transition map $f_{q}$ is unique. The RDS of the embedded
Markov chain $(X_{n})_{n\geq 0}$ is given by the tuple $(\theta,\varphi)$ with
the cocycle map
$\varphi:\mathbb{N}_{0}\times\mathcal{Q}_{+}\times\mathbb{X}\to\mathbb{X}$
defined by

$\varphi^{n}_{q}(x)=\begin{cases}f_{\theta^{n-1}q}\circ\cdots\circ
f_{q}(x)&\textup{if }n\geq 1,\\\ x&\textup{if }n=0.\end{cases}$ (11)

Noting that for any $n\in\mathbb{N}_{0}$ we have that
$\theta^{n}q=(\theta\circ...\circ\theta)q=(q_{n},q_{n+1},...)$ holds for the
shift map $\theta$ given in (9), it is straightforward to verify that the
cocycle property (1) holds.

For each initial state $x\in\mathbb{X}$ and each $q\in\mathcal{Q}_{+}$ we
obtain the orbit of states
$(x_{n})_{n\in\mathbb{N}_{0}}=(\varphi_{q}^{n}(x))_{n\in\mathbb{N}_{0}}$ from
the random difference equation

$x_{n+1}=f_{\theta^{n}q}(x_{n}),\>n\in\mathbb{N}_{0}\quad\text{and}\quad
x_{0}=x\in\mathbb{X}.$ (12)

Given an initial state $x$, we have the relation

$\mathbb{P}(X_{n}\in A|X_{0}=x)=\mathbb{P}(\varphi^{n}_{q}(x)\in A)$

for any $A\in\mathcal{B}(\mathbb{X})$. Note that by virtue of a fixed order of
reaction indices assumed for (5) and of the explicit recursion formula (3),
the Markov chain defines our RDS in a unique way; this is generally not the
case as there may be different versions of a Markov chain purely characterized
by its transition probabilities; for a general discussion see also [18].

###  2.3 RDS formulation of the augmented Markov chain

The RDS $(\theta,\varphi)$ introduced before captures only the states $X_{n}$
entered by the reaction jump process at the jump times $T_{n}$. In the
following, we formulate another RDS which takes account also of the jump times
$T_{n}$ by considering the augmented Markov chain, see (2)-(3).

The state space of the augmented Markov chain is given by
$\mathbb{X}\times[0,\infty)$ with the $\sigma$-algebra given by
$\mathcal{P}(\mathbb{X})\otimes\mathcal{B}([0,\infty))$, where
$\mathcal{P}(\mathbb{X})$ denotes the power set of $\mathbb{X}$. In analogy to
$f_{q}$ defined in (10), we consider for any $r\in\mathcal{Q}_{+}$ the mapping
$g_{r}:\mathbb{X}\times[0,\infty)\to[0,\infty)$ with

$g_{r}(x,t):=t+\tau(x,r_{0})$ (13)

for $\tau$ given in (4), such that $T_{n+1}=g_{\theta^{n}r}(X_{n},T_{n})$ in
(2). However, in contrast to $f_{q}$, this mapping depends not only on
$t\in[0,\infty)$ but also on the state $x\in\mathbb{X}$. So, the corresponding
cocycle map has to depend on state and time, as well as on
$r\in\mathcal{Q}_{+}$ and $q\in\mathcal{Q}_{+}$. Therefore, we introduce the
product noise space

$\Omega_{+}=\mathcal{Q}_{+}\times\mathcal{Q}_{+}=\\{\omega=(\omega_{n})_{n\in\mathbb{N}_{0}}:\omega_{n}=(q_{n},r_{n}),q_{n},r_{n}\in[0,1]\\}$

endowed with the product $\sigma$-algebra
$\sigma(\mathcal{Q}_{+})\otimes\sigma(\mathcal{Q}_{+})$ and the product
measure
$\mathbb{P}_{\Omega_{+}}=\lambda^{\mathbb{N}_{0}}\otimes\lambda^{\mathbb{N}_{0}}$.
By abusing the notation the corresponding shift map $\theta$ acts on both
entries of a $\omega\in{\Omega_{+}}$:

$\theta\omega=\theta(q_{n},r_{n})_{n\in\mathbb{N}_{0}}=(q_{n+1},r_{n+1})_{n\in\mathbb{N}_{0}}.$

The transformation/time-one mapping for the augmented Markov chain is given by
$h_{\omega}:\mathbb{X}\times[0,\infty)\to\mathbb{X}\times[0,\infty)$ with

$h_{\omega}(x,t):=(f_{q}(x),g_{r}(x,t)),$ (14)

where $f_{q}$ and $g_{r}$ are given in (10) and (13). The cocycle map
$\psi:\mathbb{N}_{0}\times\Omega\times\mathbb{X}\times[0,\infty)\to\mathbb{X}\times[0,\infty)$
is given by

$\psi^{n}_{\omega}(x,t)=\begin{cases}h_{\theta^{n-1}\omega}\circ\cdots\circ
h_{\omega}(x,t)&\textup{if }n\geq 1,\\\ (x,t)&\textup{if }n=0,\end{cases}$

and fulfills the cocycle property (1). We obtain

$\mathbb{P}_{\Omega_{+}}\Big{(}(X_{n},T_{n})\in
A\;\big{|}\;(X_{0},T_{0})=(x,t)\Big{)}=\mathbb{P}_{\Omega_{+}}(\psi_{\omega}^{n}(x,t)\in
A)$

for $A\in\mathcal{P}(\mathbb{X})\otimes\mathcal{B}\big{(}[0,\infty)\big{)}$
and a given initial state $x$ and starting time $t$.

We note that the first component of $\psi_{\omega}^{n}$ coincides with the
cocycle map $\varphi_{q}^{n}$ of the embedded Markov chain, i.e. we have
$(\psi_{\omega}^{n}(x,t))_{1}=\varphi_{q}^{n}(x)$ for $\varphi_{q}^{n}$ given
in (11) and $q_{n}=(\omega_{n})_{1}$, while the second component
$(\psi_{\omega}^{n}(x,t))_{2}$ referring to the time points cannot be
considered separately as a cocycle.

##### Continuous-time process realizations

By means of the RDS $\psi_{\omega}^{n}$ of the augmented Markov chain, we can
introduce a version of the continuous-time Markov jump process $(X(t))_{t\geq
0}$ starting at time $t_{0}=0$ in $X(0)=x_{0}$ by

$\Phi^{t}_{\omega}(x_{0}):=\varphi^{n}_{q}(x_{0})\quad\mbox{for}\quad(\psi^{n}_{\omega}(x_{0},0))_{2}\leq
t<(\psi^{n+1}_{\omega}(x_{0},0))_{2}.$ (15)

for $\omega=(q,r)$. This is helpful for illustrating the dynamics: In Fig 2,
common noise realizations of the continuous-time birth-death process given in
Example 1 are depicted for different initial states $x_{0}\neq y_{0}$. As we
can see, the two realizations seem to approach each other – with a certain
time-delay – given that the difference $x_{0}-y_{0}$ of the initial states is
even (see Fig 2(a)), while this is not the case for an odd difference
$x_{0}-y_{0}$ (see Fig 2(b)). This observation motivates to formulate and
analyze the synchronization behavior of random dynamical systems for the
reaction systems under consideration, which we will do in the following
section.

Importantly, we note that $\Phi^{t}_{\omega}(x_{0})$, as given in (15), does
not satisfy the cocycle property and, hence, the continuous-time Markov jump
process is itself not an RDS in this formulation. The easiest way to observe
this is that there are instances of
$\Phi^{t}_{\omega}(x_{0})=\Phi^{t}_{\omega}(y_{0})$ but
$\Phi^{t+s}_{\omega}(x_{0})\neq\Phi^{t+s}_{\omega}(y_{0})$ for some $t,s>0$,
$x_{0}\neq y_{0}$. However, the information from our RDS analysis of the
augmented Markov chain is insightful in terms of characterizing the
continuous-time process, as illustrated in Figure 2 and the following results
on time-shifted synchronization. We additionally emphasize that our
construction illustrates an intriguing lack of commutativity in the following
sense: the Markov jump process admits a version that corresponds to the
augmented Markov chain which directly induces an RDS. This RDS can be related
back to the original process via (15) giving a version of the Markov jump
process which, however, does not satisfy the cocycle property and is therefore
not part of a continuous-time RDS itself. In summary, the RDS structure lies
in the space-time version of the reaction rate process, revealing also
relevant information about this process as we will see in the following.

(a)

(b)

Figure 2: Continuous-time realizations $\Phi^{t}_{\omega}(x_{0})$ given in
(15) for the birth-death process of Example 1, driven by the same noise
$\omega$. Realizations for two initial values with (a) even distance and (b)
odd distance. In (a), the orange trajectory seems to become a time-delayed
copy of the blue one, while this is not the case in (b). The rate constants
are chosen as $\gamma_{1}=10$, $\gamma_{2}=1$.

##  3 Synchronization of reaction jump processes

In the following, we introduce the terms synchronization and partial
synchronization for the random dynamical systems under consideration. We
analyze the synchronizing properties of the birth-death process given in
Example 1 as well as of the Schlögl model of Example 2.

###  3.1 General formulation

Analogously to [16], we say that an RDS $(\theta,\varphi)$ on
$\mathbb{X}=\mathbb{N}_{0}^{L}$ is _synchronizing in $S\subset\mathbb{X}$_
(or, simply, _synchronizing_ when $S=\mathbb{X}$) if for every two different
initial states $x,y\in S$ and $\mathbb{P}$-a.e. $q\in\mathcal{Q}_{+}$ there
exists a number $n_{0}\equiv n_{0}(x,y,q)\in\mathbb{N}_{0}$ such that

$\varphi_{q}^{n_{0}}(x)=\varphi_{q}^{n_{0}}(y).$ (16)

It follows from the cocycle property that if (16) holds for some
$n_{0}\in\mathbb{N}$, it is true for any other $n\geq n_{0}$. We say that the
RDS $(\theta,\varphi)$ is partially synchronizing if there exists a partition
$\xi=\\{W_{0},\dotsc,W_{p-1}\\}$ of $\mathbb{X}$ such that $(\theta,\varphi)$
is synchronizing in each $W_{i}\in\xi$.

For a fixed $q\in\mathcal{Q_{+}}$ and two different initial states
$x,y\in\mathbb{X}$ the process
$(\varphi_{q}^{n}(x),\varphi_{q}^{n}(y))_{n\in\mathbb{N}_{0}}$ in the product
space $\mathbb{X}^{2}$ is called the two-point motion. Let $\Delta$ denote the
diagonal in $\mathbb{X}^{2}$, i.e.

$\Delta:=\\{(x,y)\in\mathbb{X}^{2}:x=y\\}.$ (17)

Hence the RDS is synchronizing if and only if the two-point motion reaches the
diagonal $\Delta$ at a time index $n_{0}(x,y,q)$.

######  Remark 1.

Since $\mathbb{X}$ is discrete, in order to show that an RDS
$(\theta,\varphi)$ is synchronizing in $S$ it suffices to show that for every
two initial states $x,y\in S$ there is a full probability measurable set
$\mathcal{Q}_{x,y}$ in which (16) holds. Indeed, for each $x,y\in\mathbb{X}$
consider the measurable sets

$\mathcal{Q}_{x,y}=\bigcup_{n=1}^{\infty}\\{q\in\mathcal{Q}_{+}:\varphi^{n}_{q}(x)=\varphi^{n}_{q}(y)\\}.$

and assume that $\mathbb{P}(\mathcal{Q}_{x,y})=1$. We can thus take
$\mathcal{Q}_{S}=\bigcap_{x,y\in S}\mathcal{Q}_{x,y}$ and (16) holds for every
$x,y\in\mathbb{X}$ and $q\in\mathcal{Q}_{S}$, where
$\mathbb{P}(\mathcal{Q}_{S})=1$.

Furthermore, if an RDS is partially synchronizing then for each $W_{i}\in\xi$
we consider the corresponding sets $\mathcal{Q}_{W_{i}}$. By considering
$\hat{Q}=\bigcap_{i=0}^{p-1}Q_{W_{i}}$, we can always assume without loss of
generality that the set is the same for each element of the partition.

##### Time-shifted synchronization

We observe that the synchronization of the RDS $(\theta,\varphi)$ of the
embedded Markov chain directly implies a time-shifted synchronization of the
RDS $(\theta,\psi)$ of the augmented Markov chain in the following sense. For
every two different initial states $x,y\in S$, an initial time
$t\in[0,\infty)$ and $\mathbb{P}_{\Omega_{+}}$-a.e. $\omega\in\Omega$ there
exists a number $n_{0}:=n_{0}(x,y,t,\omega)\in\mathbb{N}_{0}$ and a value
$R:=R(x,y,t,\omega)\geq 0$ such that

$\begin{cases}\varphi_{q}^{n}(x)=\varphi_{q}^{n}(y)\\\
\big{|}(\psi_{\omega}^{n}(x,t))_{2}-(\psi_{\omega}^{n}(y,t))_{2}\big{|}=R\end{cases}\text{for
all }n\geq n_{0}.$ (18)

This means that from a certain time point, the states $X_{n}$ of two
realizations of the augmented Markov chain actually coincide, while for the
jump times only the differences $T_{n+1}-T_{n}$ become the same. Note that we
are considering two different initial states $x,y\in S\subset\mathbb{X}$, but
start with both at the same initial time $t\in[0,\infty)$. Usually this
initial time is $t=0$. If (18) holds for all $x,y$ in each component of a
partition of $\mathbb{X}$ we analogously speak of partial time-shifted
synchronization.

The fact that the (partial) time-shifted synchronization of the RDS
$(\theta,\psi)$ follows from the (partial) synchronization of the RDS
$(\theta,\varphi)$ of the related embedded Markov chain is because the first
component $f_{q}(x)$ of the time-one mapping $h_{\omega}$ is independent of
$t$, see (14). As a consequence, for analyzing the synchronization properties
of the augmented Markov chain (and with it the synchronization properties of
the continuous-time reaction jump process) it suffices to consider the
corresponding embedded Markov chain.

We proceed by analyzing the synchronization properties for the special case of
the birth-death process.

###  3.2 Synchronization of the birth-death process

We consider the embedded Markov chain of the birth-death process defined in
Example 1, which for simplicity we refer to it simply as birth-death chain in
the following333One should bear in mind that in many references (see e.g. [8])
a birth-death chain is a general Markov chain in $\mathbb{N}_{0}$ where the
only possible transitions from $x\in\mathbb{N}$ are $x+1$ or $x-1$. In this
sense the embedded Markov chain of the Schlögl model is also a birth-death
chain but we will distinguish the two by this choice of terminology.. Its
transition map $f_{q}:\mathbb{N}_{0}\to\mathbb{N}_{0}$ can be chosen as

$f_{q}(x)=\begin{cases}x+1\quad\text{if
}q_{0}<\frac{\gamma_{1}}{\gamma_{1}+\gamma_{2}x},\\\
x-1\quad\text{otherwise,}\end{cases}$

for $q\in\mathcal{Q}_{+}$, which directly follows from Eq. (6). Recall that
the choice of $f_{q}$ is not unique as long as the order of reaction indices
it not fixed in (5). A reordering of the reactions would lead to
$f_{q}(x)=x-1$ if $q_{0}<\frac{\gamma_{2}x}{\gamma_{1}+\gamma_{2}x}$ and
$f_{q}(x)=x+1$ otherwise, defining an RDS with the same statistics and, by
symmetry, also the same topological properties that we will study throughout
the rest of the paper.

Our goal is to show that the RDS of the birth-death chain partially
synchronizes. For this purpose, we consider the two-point motion
$(\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0}))_{n\in\mathbb{N}_{0}}$ of the
birth-death chain, which is depicted in Figure 3 for two different pairs of
initial states. As a first step, we will make clear through the transition
probabilities of the two-point motion that the thick diagonal

$\mathbb{D}:=\big{\\{}(x,y)\in\mathbb{N}_{0}^{2}:y\in\\{x-1,x,x+1\\}\big{\\}}$
(19)

is forward invariant for the two-point motion, i.e.
$(\varphi^{n}_{q}(x),\varphi^{n}_{q}(y))\in\mathbb{D}$ for all
$n\in\mathbb{N}$ if $(x,y)\in\mathbb{D}$; see Figure 4 for an illustration. We
prove partial synchronization for the birth-death chain via Lemma 1,
Proposition 2, and Corollary 3.

(a)

(b)

Figure 3: Two-point motion
$(\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0}))_{n\in\mathbb{N}_{0}}$ of the
birth-death process. The trajectory in (a) belongs to the realization shown in
Figure 2(a), while (b) refers to Figure 2(b). The rate constants are chosen as
$\gamma_{1}=10$, $\gamma_{2}=1$.

##### Transition probabilities

Let $z$ be a variable which takes the values $\nu_{1}=1$ or $\nu_{2}=-1$ of
the corresponding state-change vectors. For the transition probabilities of
the RDS of the embedded Markov chain we set

$P_{z}(x):=\mathbb{P}\Big{(}\varphi_{q}^{n+1}(x_{0})=x+z\;\Big{|}\;\varphi_{q}^{n}(x_{0})=x\Big{)}$
(20)

for an arbitrary state $x\in\mathbb{N}_{0}$ and an initial state
$x_{0}\in\mathbb{N}$. This probability is independent of $n$ because the
process is time-homogeneous. Using again Eq. (6), we obtain

$\displaystyle P_{1}(x)=\frac{\gamma_{1}}{\gamma_{1}+\gamma_{2}x},\quad
P_{-1}(x)=\frac{\gamma_{2}x}{\gamma_{1}+\gamma_{2}x}.$ (21)

Just as for the one-point motion
$(\varphi_{q}^{n}(x_{0}))_{n\in\mathbb{N}_{0}}$, we can also determine the
transition probabilities for the Markovian dynamics of the two-point motion
$(\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0}))_{n\in\mathbb{N}_{0}}$. For
$z_{1},z_{2}\in\\{1,-1\\}$ set

$\displaystyle P_{(z_{1},z_{2})}(x,y):=\mathbb{P}\Big{(}\big{(}$
$\displaystyle\varphi_{q}^{n+1}(x_{0}),\varphi_{q}^{n+1}(y_{0})\big{)}=(x+z_{1},y+z_{2})$
$\displaystyle\Big{|}\big{(}\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0})\big{)}=(x,y)\Big{)}$

for any $(x,y)\in\mathbb{N}_{0}^{2}$. Given (21) we can deduce that the
transition probabilities are

$\displaystyle\hskip 56.9055ptP_{(1,1)}(x,y)=\min\\{P_{1}(x),P_{1}(y)\\},$
(yellow) (22) $\displaystyle\hskip
56.9055ptP_{(-1,1)}(x,y)=\max\\{0,P_{1}(y)-P_{1}(x)\\},$ (red) (23)
$\displaystyle\hskip 56.9055ptP_{(1,-1)}(x,y)=\max\\{0,P_{1}(x)-P_{1}(y)\\},$
(blue) (24) $\displaystyle\hskip
56.9055ptP_{(-1,-1)}(x,y)=1-\max\\{P_{1}(x),P_{1}(y)\\}.$ (green) (25)

The colors refer to the transitions given by the arrows in Figure 4. Since we
have $P_{(-1,1)}(x,y)=0$ for $x<y$ and $P_{(1,-1)}(x,y)=0$ for $x>y$, it
follows inmediately that the set $\mathbb{D}$ is forward-invariant for the
two-point motion.

##### First hitting time of $\mathbb{D}$

We will show that the absorbing set $\mathbb{D}$ is reached by the two-point
motion $(\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0}))_{n\in\mathbb{N}_{0}}$
of the birth-death process almost surely in finite time regardless of the
starting point. Formally speaking, by considering

$\tau_{\mathbb{D}}(x_{0},y_{0},q):=\inf\\{n\geq
0:(\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0}))\in\mathbb{D}\\}$

as the first hitting time of the thick diagonal $\mathbb{D}$, we show that

$\mathbb{P}\big{(}\tau_{\mathbb{D}}(x_{0},y_{0},q)<\infty\big{)}=1$

holds for all $x_{0},y_{0}\in\mathbb{N}_{0}$. To do so, we first define for a
given $d\in\mathbb{Z}$ the level set

$I_{d}:=\left\\{(x,y)\in\mathbb{N}_{0}^{2}:x-y=d\right\\}$ (26)

and show in the following lemma that for $d\neq 0$ a process starting in
$I_{d}$ will almost surely leave this set in finite time.

Figure 4: Schematic illustration for the transitions of the two-point motion
for the birth-death process. The thick diagonal $\mathbb{D}$ is defined in
(19), the level sets $I_{d}$ are defined in (26). The arrows indicate the
directions in which the two-point motion can move, with the colors indicating
the values of the corresponding probabilities as given in (22)-(25).

######  Lemma 1.

Let $d\in\mathbb{Z}\setminus\\{0\\}$ be given. Then, for each initial state
$(x_{0},y_{0})\in I_{d}$ the two-point motion
$(\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0}))_{n\in\mathbb{N}_{0}}$ of the
birth-death chain exits $I_{d}$ $\mathbb{P}$-a.s., i.e.,

$\mathbb{P}\big{(}(\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0}))\in
I_{d}\;\forall n\geq 0\big{)}=0\quad\text{for all }(x_{0},y_{0})\in I_{d}.$

###### Proof.

Without loss of generality we only consider the case in which $d\geq 1$. For
$x\in\mathbb{N}_{0}$ consider the state $(x+d,x)\in I_{d}$ and define

$p_{x}:=\mathbb{P}\big{(}(\varphi_{q}^{n}(x+d),\varphi_{q}^{n}(x))\in
I_{d}\;\forall n\geq 0\big{)}$

as the probability for the two-point motion to stay forever on $I_{d}$ given
that it starts in $(x+d,x)\in I_{d}$. By means of the law of total probability
we have

$p_{x}=P_{(1,1)}(x+d,x)\cdot p_{x+1}+P_{(-1,-1)}(x+d,x)\cdot p_{x-1}$

for $x\geq 1$, where $P_{(1,1)}(x+d,x)=\frac{1}{1+\alpha(x+d)}$ and
$P_{(-1,-1)}(x+d,x)=\frac{\alpha x}{1+\alpha x}$ with
$\alpha:=\frac{\gamma_{2}}{\gamma_{1}}$, see (22) and (25). From this we can
deduce the second-order difference equation

$p_{x+2}=(1+\alpha(x+d+1))\left(p_{x+1}-\frac{\alpha(x+1)}{1+\alpha(x+1)}p_{x}\right)$
(27)

for $x\geq 0$. Moreover, for $x=0$ we have $p_{0}=P_{(1,1)}(d,0)\cdot p_{1}$
such that

$p_{1}=(1+\alpha d)p_{0}.$ (28)

Since $p_{1}$ is proportional to $p_{0}$, it follows inductively from Eq. (27)
that $p_{x}$ is proportional to $p_{0}$ for all $x\in\mathbb{N}_{0}$. We prove
by contradiction that the sequence $(p_{x})_{x\in\mathbb{N}_{0}}$ of
probabilities has to fulfill $p_{x}=0$ for all $x\in\mathbb{N}_{0}$.

Indeed, assume that $p_{0}>0$. From (28) we obtain

$p_{1}-p_{0}=\left(1-\frac{1}{1+\alpha d}\right)p_{1}=\alpha dp_{0}>0.$ (29)

On the other hand, it follows immediately from (27) that

$p_{x+2}-p_{x+1}=\alpha(x+d+1)p_{x+1}-\frac{\alpha(x+1)\left[1+\alpha(x+d+1)\right]}{1+\alpha(x+1)}p_{x}.$

Hence, by adding and subtracting
$\frac{\alpha(x+1)\left[1+\alpha(x+d+1)\right]}{1+\alpha(x+1)}p_{x+1}$ we
obtain

$\displaystyle p_{x+2}-p_{x+1}$ $\displaystyle=\frac{\alpha
d}{1+\alpha(x+1)}\cdot
p_{x+1}+\frac{\alpha(x+1)\left[1+\alpha(x+d+1)\right]}{1+\alpha(x+1)}\cdot(p_{x+1}-p_{x})$
$\displaystyle\geq\alpha(x+1)\left(1+\frac{\alpha
d}{1+\alpha(x+1)}\right)(p_{x+1}-p_{x}).$

Let $u_{x}:=p_{x+1}-p_{x}$. Then, from the last inequality it follows that

$u_{x+1}>\alpha(x+1)u_{x}.$

By iterating the above inequality it follows that

$u_{x}>\alpha^{x-1}\cdot x!\;u_{0}\qquad x\in\mathbb{N}.$

Since $u_{0}>0$, see (29), and since $\alpha^{x-1}\cdot x!\rightarrow\infty$
as $x\rightarrow\infty$ for any $\alpha>0$, this implies that
$u_{x}\rightarrow\infty$. It then follows that $p_{x}\rightarrow\infty$, which
is a contradiction since $p_{x}\in[0,1]$. In conclusion, we obtain $p_{0}=0$,
and by proportionality, $p_{x}=0$ for all $x$. Noticing that all states in
$I_{d}$ are of the form $(x+d,x)$ for some $x\in\mathbb{N}_{0}$ completes the
proof. ∎

By means of Lemma 1 we can now make the following central statement.

######  Proposition 2.

For each $(x,y)\in\mathbb{N}_{0}^{2}$, the two-point motion
$(\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0}))_{n\in\mathbb{N}_{0}}$ of the
birth-death chain reaches the thick diagonal $\mathbb{D}$ almost surely in
finite time.

###### Proof.

Let $(x_{0},y_{0})\in I_{d}$ for a given $d=x_{0}-y_{0}\neq 0$. According to
Lemma 1, the two-point motion
$(\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0}))_{n\in\mathbb{N}_{0}}$ almost
surely escapes from $I_{d}$ in finite time. Given the transition probabilities
(22)-(25), it can only end up in $I_{d-2}$ when $d\geq 1$ or in $I_{d+2}$ when
$d\leq-1$. This happens a finite number of times until the process reaches
$\mathbb{D}=I_{-1}\cup I_{0}\cup I_{1}$. ∎

######  Corollary 3.

The RDS $(\theta,\varphi)$ for the embedded Markov chain of the birth-death
process partially synchronizes, and the partition is given by
$\xi=\\{W_{0},W_{1}\\}$ with $W_{0}=\\{0,2,4,\ldots\\}$ and
$W_{1}=\\{1,3,5,\ldots\\}$.

###### Proof.

Let $W_{0},W_{1}$ be as in the statement above. As observed before in
Proposition 2, for any fixed $(x,y)\in I_{d}$ with $d\neq 0$ the two-point
motion escapes in finite time to $I_{d-2}$ if $d\geq 1$, or to $I_{d+2}$ when
$d\leq-1$. Thus, it reaches $\Delta=I_{0}$ in finite time $\mathbb{P}$-a.s. if
and only if $d$ is even, that is if $x,y\in W_{i}$ for some $i\in\\{0,1\\}$.
The result follows from Proposition 2 and Remark 1. ∎

Corollary 3 states that whenever $x,y\in W_{i}$ ($i=0,1$) we have for almost
all $q\in\mathcal{Q}_{+}$ that $\\#\varphi^{n}_{q}(\\{x,y\\})=1$ for all $n$
sufficiently large, where $\\#A$ denotes the cardinality of a set $A$. More
generally, for each finite (deterministic) set $K\subset W_{i}$ we obtain that
for almost all $q\in\mathcal{Q}_{+}$ there is a $n_{0}(K,q)\in\mathbb{N}$ such
that $\\#\varphi^{n}_{q}(K)=1$ for all $n\geq n_{0}(K,q)$. This almost sure
convergence implies the convergence in probability given by

$\lim_{n\rightarrow\infty}\mathbb{P}(\\#\varphi^{n}_{q}(K)\geq 2)=0.$ (30)

This last statement can be extended to finite random sets as defined in the
following Definition 4. For this, let
$d:\mathbb{X}\times\mathbb{X}\to[0,\infty)$ be the Euclidean distance on
$\mathbb{X}$ and define

$d(x,B):=\inf_{y\in B}d(x,y)$ (31)

for non-empty sets $B\subset\mathbb{X}$.

######  Definition 4.

Let $(\mathcal{Q},\sigma(\mathcal{Q}),\mathbb{P})$ be an arbitrary probability
space and $\mathbb{X}=\mathbb{N}_{0}^{L}$. A mapping
$K:\mathcal{Q}\rightarrow\mathcal{P}(\mathbb{X})$, denoted as $q\mapsto
K_{q}$, is a random set if the function $q\mapsto d(x,K_{q})$ is measurable
for each $x\in\mathbb{X}$.

We say that a random set $K:\mathcal{Q}\rightarrow\mathcal{P}(\mathbb{X})$ is
a finite random set if $K_{q}$ is nonempty and finite for every
$q\in\mathcal{Q}$. A finite random set in $\mathcal{P}(\mathbb{X})$ is
contained in a deterministic finite set with high probability as indicated in
the next proposition.

######  Proposition 5.

Let $K:\mathcal{Q}\rightarrow\mathcal{P}(\mathbb{X})$ be a finite random set.
Then, for each $\varepsilon>0$ there is a finite set
$F_{\varepsilon}\subset\mathcal{P}(\mathbb{X})$ such that

$\mathbb{P}\left(K_{q}\subset F_{\varepsilon}\right)\geq 1-\varepsilon.$

###### Proof.

The statement is a particular case of a more general setting, see [5,
Proposition 3.15]. ∎

We can now generalize property (30) for the birth-death chain for arbitrary
random finite sets in the next Proposition.

Figure 5: Time-shifted synchronization for the birth-death process. Extract of
Figure 2(a), showing two realizations of $\Phi_{\omega}^{t}(x_{0})$ for
initial values with even distance. The time-shifted synchronization starts at
time index $n_{0}=19$, where both realizations reach the state $x=8$. For the
process starting in $x_{0}=5$ this happens at time
$T_{\text{syn}}=(\psi_{\omega}^{n_{0}}(5,0))_{2}\approx 1.395$, while for the
other realization the time point is given by
$T^{\prime}_{\text{syn}}=(\psi_{\omega}^{n_{0}}(15,0))_{2}\approx 1.064$.
I.e., the time-shift is $R=T_{\text{syn}}-T^{\prime}_{\text{syn}}\approx
0.331$ for these initial states.

######  Proposition 6.

Consider the setting of the birth-death chain. Let
$K:\mathcal{Q}_{+}\rightarrow\mathcal{P}(\mathbb{N}_{0})$ be a random finite
set such that $K_{q}\subset W_{i}$ $\mathbb{P}$-a.s. for some $i\in\\{0,1\\}$.
Then $\\#\varphi^{n}_{q}(K_{q})\rightarrow 1$ in probability.

###### Proof.

Since $K_{q}$ is nonempty $\mathbb{P}$-a.s., the set $\varphi^{n}_{q}(K_{q})$
has at least one element for all $n\in\mathbb{N}$ $\mathbb{P}$-a.s.. On the
other hand, let $\varepsilon>0$ be arbitrarily small and
$F_{\varepsilon}\subset\mathbb{N}_{0}$ a finite set as in Proposition 5. Then,

$\mathbb{P}(\\#\varphi^{n}_{q}(K_{q})\geq
2)=\mathbb{P}\left(\\#\varphi^{n}_{q}(K_{q})\geq 2\ \land K_{q}\subset
F_{\varepsilon}\right)+\mathbb{P}\left(\\#\varphi^{n}_{q}(K_{q})\geq 2\ \land
K_{q}\not\subset F_{\varepsilon}\right).$

Note that for any $n\in\mathbb{N}_{0}$,

$\\{q\in\mathcal{Q}_{+}:\\#\varphi^{n}_{q}(K_{q})\geq 2\ \land\ K_{q}\subset
F_{\varepsilon}\\}\subset\\{q\in\mathcal{Q}_{+}:\\#\varphi^{n}_{q}(F_{\varepsilon})\geq
2\\},$

and thus we deduce together with Proposition 5 that

$\mathbb{P}(\\#\varphi^{n}_{q}(K_{q})\geq
2)\leq\mathbb{P}(\\#\varphi^{n}_{q}(F_{\varepsilon})\geq 2)+\varepsilon.$

Due to (30) and since $\varepsilon$ was arbitrarily small we conclude that

$\lim_{n\rightarrow\infty}\mathbb{P}(\\#\varphi^{n}_{q}(K_{q})\geq 2)=0,$

and the result follows. ∎

As noted in the end of Sec. 3.1, the partial synchronization of the RDS
$(\theta,\varphi)$ for the embedded Markov chain directly implies the partial
time-shifted synchronization of the corresponding augmented Markov process,
which can easily be seen from Eq. (2). Thus, the observations from Figure 2
can now be confirmed/clarified: In Figure 2(a), we have $x_{0}\in W_{0}$ for
both initial states, such that time-shifted synchronization as defined in (18)
is guarantied by Corollary 3, see Figure 5 for a detailed look at the
dynamics. In contrast, the initial states chosen in Figure 2(b) are not in the
same set $W_{i}$, and consequently, the trajectories do not synchronize.
However, they are likely to stay close to each other because the corresponding
two-point motion reaches the thick diagonal where the distance between states
is not larger than one.

###  3.3 Synchronization for the Schlögl model

In this section, we consider Example 2 as a variation of the birth-death
process with additional bistable structure. We analyze the RDS
$(\theta,\varphi)$ for the embedded Markov chain of the according reaction
jump process. With the notation as in (20), the transition probabilities of
the system are

$P_{1}(x)=\frac{\gamma_{1}+\gamma_{3}x(x-1)}{\mu(x)},\quad
P_{-1}(x)=\frac{\gamma_{2}x+\gamma_{4}x(x-1)(x-2)}{\mu(x)}$ (32)

for an arbitrary state $x\geq 0$. In Figure 6, two continuous-time
realizations of the dynamics are shown. As in the birth-death process, one can
observe a time-shifted synchronization when starting with even distance, see
Figure 6(a). For an odd distance in the starting points, see Figure 6(b), the
separation of the trajectories is even more significant than in the birth-
death scenario.

(a)

(b)

Figure 6: Continuous-time realizations $\Phi^{t}_{\omega}(x_{0})$ given in
(15) for the Schlögl model of Example 2, driven by the same noise.
Realizations for two initial values with (a) even distance and (b) odd
distance. In (a), the orange trajectory becomes a time-delayed copy of the
blue one, just as in Figure 2, while in (b), the trajectories clearly
separate. The rate constants are chosen as $\gamma_{1}=6$,
$\gamma_{2}=3.5,\gamma_{3}=0.4,\gamma_{4}=0.0105$.

The transition probabilities of the corresponding two-point motion are
analogous to the ones for the birth-death process given in (22) - (25),
however with $P_{1}$ and $P_{-1}$ defined in (32). Figure 7 shows two
realizations of the two-point motion for different initial states. As opposed
to the birth-death chain, the thick diagonal $\mathbb{D}$ is not absorbing in
the Schlögl model. This can be seen from the transition probabilities of the
two-point motion which are non-zero also for directions pointing away from the
thick diagonal (see Figure 8).

(a)

(b)

Figure 7: Two-point motion
$(\varphi_{q}^{n}(x_{0}),\varphi_{q}^{n}(y_{0}))_{n\in\mathbb{N}_{0}}$ of the
Schlögl model. The trajectory in (a) belongs to the realization shown in
Figure 6(a), while (b) refers to Figure 6(b).

A proof of the conjectured partial synchronization is much more involved than
in the simple birth-death case, due to the described lack of monotonicity in
the two-point dynamics towards the diagonal. Hence, the proof strategy of
Proposition 2 and the following corollaries does not apply here. We leave a
suitably generalized proof with different means as an open problem for the
future. From such a partial synchronization of the embedded Markov chain, one
may conclude directly the time-shifted synchronization of the continuous-time
realization $(\Phi^{t}_{\omega}(x))_{t\geq 0}$ for initial states with even
distance as observed in Figure 6(a).

Figure 8: Possible transitions to the top left (red) and the bottom right
(blue) for the Schlögl model. Transitions parallel to the diagonal are not
displayed.

##  4 Random attractors associated to (embedded) Markov chains

In this section we introduce different notions of random attractors for an
RDS, providing insights into their relationships in the context of a discrete
state in Sec. 4.1. As a matter of fact, we show in Sec. 4.2 that under very
mild conditions an RDS induced by the random difference equation (12) admits a
weak attractor. The partial synchronization for the embedded Markov chain of
the birth-death process, which was proved in Section 3, will be useful to
describe in full detail the structure of its attractor and the related sample
measures in Sec. 4.3. We explore numerically such characteristics for the
attractor of the embedded Markov chain of the Schlögl model in Sec. 4.4.

###  4.1 General properties of random attractors in discrete time and discrete
space

In this and the following section (Sec. 4.2) we consider a generalized setting
as compared to Sec. 2.2, since our insights on random attractors for the
discrete case can be easily formulated in that generality. More concretely,
let $\mathcal{Q}$ be a probability space with measure $\mathbb{P}$ and let
$\theta:\mathcal{Q}\to\mathcal{Q}$ be an invertible and $\mathbb{P}$-ergodic
map. Furthermore, for $\mathbb{X}=\mathbb{N}_{0}^{L}$ as before, let
$\varphi:\mathbb{N}_{0}\times\mathcal{Q}\times\mathbb{X}\to\mathbb{X}$ be a
measurable cocycle (cf. property (1)) such that $(\theta,\varphi)$ is a random
dynamical system. One of the characteristics of random attractors is that of
invariance as described in the following.

######  Definition 7.

Given a RDS $(\theta,\varphi)$, we say that a random set
$A:\mathcal{Q}\rightarrow\mathcal{P}(\mathbb{X})$ (recall Definition 4) is
invariant (or $\varphi$-invariant) if

$\varphi_{q}^{n}(A_{q})=A_{\theta^{n}q}\qquad\textup{ for all
}n\in\mathbb{N},\ \mathbb{P}\textup{-a.s.}$ (33)

Notice that without loss of generality, we can assume that (33) is satisfied
everywhere by restricting ourselves to the full probability set where this is
true.

######  Remark 2.

Since we deal with a discrete-time setting, condition (33) is fulfilled as
soon as it is satisfied for $n=1$, that is $\varphi_{q}^{1}(A_{q})=A_{\theta
q}$ $\mathbb{P}$-a.s. Indeed, let $\mathcal{Q}_{1}\subset\mathcal{Q}$ be a
measurable set of full probability such that (33) holds for $n=1$. By the
invariance of $\mathbb{P}$ under $\theta$, we have that
$\mathbb{P}(\mathcal{Q}_{1})=\mathbb{P}(\theta^{-1}\mathcal{Q}_{1})=1$.
Consider then the full probability set
$\mathcal{Q}_{2}=\theta^{-1}\mathcal{Q}_{1}\cap\mathcal{Q}_{1}$, which is the
set where (33) holds for $n=1,2$, and inductively for $k\in\mathbb{N}$ the set
$\mathcal{Q}_{k+1}=\theta^{-1}\mathcal{Q}_{k}\cap\mathcal{Q}_{k}$ for which
(33) is satisfied for all $n\leq k+1$. Hence, condition (33) holds on
$\tilde{\mathcal{Q}}:=\bigcap_{n=1}^{\infty}\mathcal{Q}_{n}$, and
$\mathbb{P}(\tilde{\mathcal{Q}})=1$.

We give the notion of attraction in terms of the Hausdorff semi-distance for
non-empty sets

$\mathrm{dist}(A,B)=\sup_{x\in A}d(x,B),\quad A,B\subset\mathbb{X},$

where $d(x,B)$ is given in (31). In general, an invariant compact random set
$A$ (i.e., $A_{q}$ is compact for all $q\in\mathcal{Q}$) is called

  * (i)

a (strong) forward attractor, if for each compact set $B\subset\mathbb{X}$

$\lim_{n\to\infty}\mathrm{dist}(\varphi_{q}^{n}(B),A_{\theta^{n}q})=0\quad\mathbb{P}\text{-a.s.},$
(34)

  * (ii)

a (strong) pullback attractor if for each compact set $B\subset\mathbb{X}$

$\lim_{n\to\infty}\mathrm{dist}(\varphi_{\theta^{-n}q}^{n}(B),A_{q})=0\quad\mathbb{P}\text{-a.s.},$

  * (iii)

a weak attractor if for each compact set $B\subset\mathbb{X}$

$\lim_{n\to\infty}\mathrm{dist}(\varphi_{q}^{n}(B),A_{\theta^{n}q})=0\ \text{
in probability,}$ (35)

cf. e.g. [7]. By the discrete nature of the state space, compact sets in
$\mathbb{X}$ are simply finite sets. Moreover, for any $q\in\mathcal{Q}$ and
$B\subset\mathbb{X}$,

$\lim_{n\to\infty}\mathrm{dist}\left(\varphi^{n}_{\theta^{-n}q}(B),A_{q}\right)=0\
\text{ iff }\ \exists N\equiv N(q,B)\text{ s.t. }\forall n\geq N,\
\varphi^{n}_{\theta^{-n}q}(B)\subset A_{q}.$ (36)

In a similar fashion,
$\lim_{n\rightarrow\infty}\mathrm{dist}(\varphi^{n}_{q}(B),A_{\theta^{n}q})=0$
if and only if there exists $N\equiv N(q,B)$ such that
$\varphi^{N}_{q}(B)\subset A_{\theta^{N}q}$. Here, the invariance of $A$
guarantees that $\varphi^{n}_{q}(B)\subset A_{\theta^{n}q}$ for all $n\geq N$.

It is known that weak attractors are unique, in the sense that two weak
attractors $A$ and $\tilde{A}$ coincide $\mathbb{P}$-a.s. [11, Lemma 1.3].
Observe also that in the weak sense “forward” and “pullback” attraction are
the same since $\theta$ is measure-preserving. In other words, $A$ is a weak
attractor if and only if for each compact set $B\subset\mathbb{X}$

$\lim_{n\to\infty}\mathrm{dist}(\varphi_{\theta^{-n}q}^{n}(B),A_{q})=0\ \text{
in probability.}$

Since the state space is discrete it is straightforward to see that if the
strong pullback convergence holds for all point sets, i.e. for $B=\\{x\\}$
where $x\in\mathbb{X}$, then it also holds for any finite $B$. Using the
terminology from [7, 6], this would imply that strong point attractors are
equivalent to strong (set) attractors as defined above. Indeed, assuming that
(36) holds for all point sets $B=\\{x\\}$, where $x\in K$ for a given finite
set $K$, it also $\mathbb{P}$-a.s. holds for $B=K$ by taking
$N(q,K)=\max\\{N(q,x):x\in K\\}$. In a similar spirit, we give equivalent
conditions for an invariant random finite set to be a weak attractor in the
next theorem.

######  Theorem 8.

Let $A:\mathcal{Q}\rightarrow\mathcal{P}(\mathbb{X})$ be an invariant compact
finite set. Then the following are equivalent.

  * I.

$A$ is a weak attractor.

  * II.

$A$ is a forward attractor.

  * III.

$A$ is a forward point attractor, that is (34) is satisfied for all
$B=\\{x\\}$ with $x\in\mathbb{X}$.

  * IV.

$A$ is a _weak point attractor_ , that is (35) is satisfied for all
$B=\\{x\\}$ with $x\in\mathbb{X}$.

  * V.

$A$ weakly attracts random finite sets, that is for any finite random set
$K:\mathcal{Q}\rightarrow\mathcal{P}(\mathbb{X})$ we have that

$\lim_{n\rightarrow\infty}\mathrm{dist}(\varphi^{n}_{q}(K_{q}),A_{\theta^{n}q})=0\
\text{ in probability.}$

###### Proof.

I implies II. Assume that $A$ is a weak attractor. Take $x\in\mathbb{X}$ as an
arbitrary initial condition and consider the sets
$\mathcal{Q}_{n}^{x}:=\\{q\in\mathcal{Q}:\varphi^{n}_{q}(x)\in
A_{\theta^{n}q}\\}$. Since the state space is discrete, we have

$\mathbb{P}(\mathcal{Q}_{n}^{x})=1-\mathbb{P}\left(d\left(\varphi^{n}_{q}(x),A_{\theta^{n}q}\right)\geq
r\right)$

for any $r\in(0,1]$. From the convergence in probability to the attractor (cf.
(35)) it follows that for an arbitrary $\varepsilon>0$ there exists $N\equiv
N(x,\varepsilon)$ such that $\mathbb{P}(\mathcal{Q}_{n}^{x})\geq
1-\varepsilon$ for all $n\geq N$. Furthermore, from the $\varphi$-invariance
of $A$ we know that $\mathcal{Q}_{n}^{x}\subset\mathcal{Q}_{n+1}^{x}$ for all
$n\in\mathbb{N}$.

Let
$\mathcal{Q}^{x}:=\bigcup_{n\in\mathbb{N}}\mathcal{Q}_{n}^{x}=\\{q\in\mathcal{Q}:\lim_{n\to\infty}d\left(\varphi^{n}_{q}(x),A_{\theta^{n}q}\right)=0\\}$.
Notice that for any $n\geq N$

$\mathbb{P}(\mathcal{Q}^{x})\geq\mathbb{P}(\mathcal{Q}_{n}^{x})\geq
1-\varepsilon,$

and since $\varepsilon$ was arbitrarily small we have that
$\mathbb{P}(\mathcal{Q}^{x})=1$. Now, we set

$\mathcal{Q}_{B}:=\bigcap_{x\in
B}\mathcal{Q}^{x}=\left\\{q\in\mathcal{Q}:\lim_{n\to\infty}d\left(\varphi^{n}_{q}(x),A_{\theta^{n}q}\right)=0\;\forall
x\in B\right\\}.$

As $B$ is finite, we have

$\mathcal{Q}_{B}=\left\\{q\in\mathcal{Q}:\lim_{n\to\infty}\mathrm{dist}\left(\varphi^{n}_{q}(B),A_{\theta^{n}q}\right)=0\right\\},$

and, using $\mathbb{P}(\mathcal{Q}^{x})=1$ for all $x$, we get
$\mathbb{P}(\mathcal{Q}_{B})=1$ meaning that $A$ is a forward attractor (cf.
(34)).

II implies III. This follows directly by taking $B=\\{x\\}$.

III implies IV. This statement is true from the fact that convergence a.s.
implies convergence in probability.

IV implies I. Assume that (35) holds for $B=\\{x\\}$ for any $x\in\mathbb{X}$,
which is equivalent to

$\lim_{n\rightarrow\infty}\mathbb{P}\left(\varphi^{n}_{q}(x)\not\in
A_{\theta^{n}q}\right)=0.$

Let $K=\\{x_{1},\ldots,x_{m}\\}\subset\mathbb{X}$ for some $m\in\mathbb{N}$.
Since for any $n\in\mathbb{N}$

$\left\\{q:\mathrm{dist}(\varphi^{n}_{q}(K),A_{\theta^{n}q})\geq
1\right\\}=\bigcup_{i=1}^{m}\\{q:\varphi^{n}_{q}(x_{i})\not\in
A_{\theta^{n}q}\\},$

it follows that

$\lim_{n\rightarrow\infty}\mathbb{P}\left(\mathrm{dist}(\varphi^{n}_{q}(K),A_{\theta^{n}q})\geq
1\right)\leq\sum_{i=1}^{m}{\lim_{n\rightarrow\infty}\mathbb{P}\left(\varphi^{n}_{q}(x_{i})\not\in
A_{\theta^{n}q}\right)}=0,$

and the claim follows.

I if and only if V. Clearly, V implies I, since for each finite set
$B\subset\mathbb{X}$ one can take $K_{q}=B$ for all $q$. Conversely, if $K$ is
a finite random set and $\varepsilon>0$ an arbitrarily small constant,
consider $F_{\varepsilon}>0$ as given in Proposition 5 so that

$\mathbb{P}\left(K_{q}\subset F_{\varepsilon}\right)\geq 1-\varepsilon.$ (37)

Note that for any $n,m\in\mathbb{N}$ we have that

$\displaystyle\mathbb{P}\left(\varphi^{n}_{q}(K_{q})\not\subset
A_{\theta^{n}}q\right)$
$\displaystyle=\mathbb{P}\left(\varphi^{n}_{q}(K_{q})\not\subset
A_{\theta^{n}}q\land K_{q}\subset
F_{\varepsilon}\right)+\mathbb{P}\left(\varphi^{n}_{q}(K_{q})\not\subset
A_{\theta^{n}}q\land K_{q}\not\subset F_{\varepsilon}\right).$

Hence, it follows from (37) and the observation

$\\{\varphi^{n}_{q}(K_{q})\not\subset A_{\theta^{n}}q\ \land\ K_{q}\subset
F_{\varepsilon}\\}\subset\\{\varphi^{n}_{q}(F_{\varepsilon})\not\subset
A_{\theta^{n}}q\\},$

that for each $\varepsilon>0$ there is $m\in\mathbb{N}$ such that

$\mathbb{P}\left(\varphi^{n}_{q}(K_{q})\not\subset
A_{\theta^{n}}q\right)\leq\mathbb{P}\left(\varphi^{n}_{q}(F_{\varepsilon})\not\subset
A_{\theta^{n}}q\right)+\varepsilon.$

Since $F_{\varepsilon}$ is deterministic and finite, statement I implies
$\mathbb{P}\left(\varphi^{n}_{q}(F_{\varepsilon})\not\subset
A_{\theta^{n}}q\right)\to 0$ for $n\to\infty$. As $\varepsilon$ was
arbitrarily small, this implies
$\mathbb{P}\left(\varphi^{n}_{q}(K_{q})\not\subset
A_{\theta^{n}}q\right)\rightarrow 0$ for $n\to\infty$. ∎

######  Remark 3.

Note that, by virtue of this theorem, weak set and point attractors are the
same for this discrete-time and discrete-space setting. Hence, the weak set or
point attractor being a random point almost surely is equivalent here; in
particular, that implies that the distinction between synchronization and weak
synchronization defined via an attractor being a singleton, as done in [11],
is not necessary here. In addition, note that, again by Theorem 8, our
definition of synchronization (Sec. 3.1) conincides with the one in [11], if
an attractor exists (under extension of $\mathcal{Q}_{+}$ to $\mathcal{Q}$).

Due to their invariance, it becomes relevant to understand the dynamics within
the attractors. In particular, the attractor for the birth-death chain admits
a periodic behaviour, as we see later in Theorem 14.

######  Definition 9.

A random set $A:\mathcal{Q}\rightarrow\mathcal{P}(\mathbb{X})$ is a _random
periodic orbit_ of period $M$ for the RDS $(\theta,\varphi)$ if for
$\mathbb{P}$-a.e. $q\in\mathcal{Q}$, $A_{q}=\\{a_{0}(q),\dotsc,a_{M-1}(q)\\}$
such that

$\varphi_{q}^{1}(a_{i}(q))=a_{i+1(\bmod M)}(\theta q)\quad\text{for
}i=0,\dotsc,M-1.$

Clearly a random periodic orbit is in particular $\varphi$-invariant.
Furthermore, we say that $A:\mathcal{Q}\rightarrow\mathcal{P}(\mathbb{X})$ is
a _(weak) attracting random periodic orbit_ if it is a random periodic orbit
and a weak attractor. Note that this can be seen as a discrete-time analogue
to the continuous-time oriented definition of a random periodic solution [38]
and its generalization [10].

###  4.2 Existence of weak attractors

In this subsection we provide general conditions for an RDS as given in the
previous section to admit a weak attractor. For this purpose, we combine [7,
Thm. 10] on the existence of a weak attractor with properties of the
stationary distribution (in case of its existence) for Markov chains on
countable state spaces (cf. [8, Chapter 6]).

######  Theorem 10.

Consider an RDS $(\theta,\varphi)$ on $\mathbb{X}$, as defined in Sec. 4.1,
and assume that for every $x\in\mathbb{X}$ the set
$\varphi^{1}_{\mathcal{Q}}(x):=\\{\varphi^{1}_{q}(x):q\in\mathcal{Q}\\}\subset\mathbb{X}$
is finite, and that the Markov chain associated to $(\theta,\varphi)$ is
irreducible and recurrent. If it admits a stationary distribution, then the
RDS admits a weak attractor.

###### Proof.

In order to prove the theorem we use the following (in fact, equivalent)
criterion for the existence of a weak attractor from [7, Thm. 10]: For every
$\varepsilon>0$ there exists a compact set $C_{\varepsilon}\subset\mathbb{X}$
such that for every compact set $K\subset\mathbb{X}$ there is a
$n_{0}\in\mathbb{N}$ so that for all $n\geq n_{0}$

$\mathbb{P}(\varphi_{q}^{n}(K)\subset C_{\varepsilon})\geq 1-\varepsilon.$

Recall that a recurrent state $x\in\mathbb{X}$ has period $\ell\geq 1$ if
$\ell$ is the greatest common denominator of the set

$\left\\{n\geq 1:\mathbb{P}\left(\varphi^{n}_{\omega}(x)=x\right)>0\right\\}.$

Since the chain is irreducible, every state has the same period. Moreover, the
state space admits a cyclic decomposition given as $\mathbb{X}=W_{0}\cup
W_{1}\cup\cdots\cup W_{\ell-1}$, where for $x\in W_{i}$ and
$i\in\\{0,1,\ldots,\ell-1\\}$

$\mathbb{P}\left(\varphi^{1}_{q}(x)\in W_{j}\right)=\begin{cases}1&\text{if
}j=i+1\ (\text{mod }\ell),\\\ 0&\text{otherwise}\end{cases},$

and the $\ell$-step process $\left(\varphi^{\ell
n}_{q}(x)\right)_{n\in\mathbb{N}_{0}}$ is aperiodic and irreducible in each
$W_{i}$ (see for instance [8, Lemma 6.7.1]). Denote by $\rho$ the unique
stationary distribution of the system and let $\varepsilon$ be arbitrarily
small. Consider $z_{\varepsilon}\in\mathbb{N}$ large enough such that

$\max_{0,1,\ldots\ell-1}\frac{\rho((\mathbb{X}\setminus C_{\varepsilon})\cap
W_{i}))}{\rho(W_{i})}<\frac{\varepsilon}{2}$ (38)

for the set $C_{\varepsilon}:=\\{0,1,2,\ldots z_{\varepsilon}\\}^{L}$.

Step 1. For each $i=0,\ldots,\ell-1$, the $\ell$-step process
$(\varphi_{q}^{\ell n}(x))_{n\in\mathbb{N}_{0}}$ starting in $x\in W_{i}$
admits a unique stationary distribution $\tilde{\rho}_{i}$ supported on
$W_{i}$.

For $i=0,\ldots,\ell-1$ and $A\subset\mathbb{X}$ set

$\tilde{\rho}_{i}(A):=\frac{\rho(A\cap W_{i})}{\rho(W_{i})}.$ (39)

Then, for any $n\in\mathbb{N}$ we have

$\displaystyle\sum_{x=0}^{\infty}{\mathbb{P}(\varphi_{q}^{\ell n}(x)\in
A)\tilde{\rho}_{i}(x)}$ $\displaystyle=\frac{1}{\rho(W_{i})}\sum_{x\in
W_{i}}{\mathbb{P}(\varphi_{q}^{\ell n}(x)\in A)\rho(x)}$
$\displaystyle\stackrel{{\scriptstyle(*)}}{{=}}\frac{1}{\rho(W_{i})}\sum_{x\in
W_{i}}{\mathbb{P}(\varphi_{q}^{\ell n}(x)\in A\cap W_{i})\rho(x)}$
$\displaystyle=\frac{1}{\rho(W_{i})}\sum_{x=0}^{\infty}{\mathbb{P}(\varphi_{q}^{\ell
n}(x)\in A\cap W_{i})\rho(x)}$ $\displaystyle=\frac{\rho(A\cap
W_{i})}{\rho(W_{i})}=\tilde{\rho}_{i}(A),$

which implies that $\tilde{\rho}_{i}$ is a stationary distribution of
$(\varphi_{q}^{\ell n}(x))_{n\in\mathbb{N}_{0}}$ on $W_{i}$. In $(*)$ we used
the fact that $\mathbb{P}(\varphi_{q}^{\ell n}(x_{0})=y)=0$ $\forall x_{0}\in
W_{i},y\notin W_{i}$. As the $\ell$-step process is irreducible in each
$W_{i}$, the stationary distribution is unique (see [8, Theorem 6.5.5]).

Step 2. Let $x\in\mathbb{X}$ and $\varepsilon>0$ be given as before. Then,
there exists $n_{1}=n_{1}(x)\in\mathbb{N}$ such that
${\mathbb{P}(\varphi_{q}^{\ell n}(x)\not\in C_{\varepsilon})<\varepsilon}$ for
all $n\geq n_{1}$.

Consider $x\in W_{i}$ for any $i\in\\{0,\ldots,\ell-1\\}$. Since the
$\ell$-step process $(\varphi_{q}^{\ell n}(x))_{n\in\mathbb{N}_{0}}$ is
irreducible and admits a stationary distribution, it is then positive
recurrent on $W_{i}$. Furthermore, since it is aperiodic in $W_{i}$,
convergence $\mathbb{P}(\varphi_{q}^{\ell
n}(x)\in\cdot)\rightarrow\tilde{\rho}_{i}$ in total variation holds as
$n\to\infty$, see [8, Theorem 6.6.4]. This implies in particular the
convergence $\mathbb{P}(\varphi_{q}^{\ell n}(x)\not\in
C_{\varepsilon})\rightarrow\tilde{\rho}_{i}(\mathbb{X}\setminus
C_{\varepsilon})$. Hence, there exists $n_{1}=n_{1}(x)\in\mathbb{N}$ such that
for any $n\geq n_{1}$ we obtain

$\mathbb{P}(\varphi_{q}^{\ell n}(x)\not\in
C_{\varepsilon})<\tilde{\rho}_{i}(\mathbb{X}\setminus
C_{\varepsilon})+\frac{\varepsilon}{2}\stackrel{{\scriptstyle\eqref{tilde_rho}}}{{=}}\frac{\rho((\mathbb{X}\setminus
C_{\varepsilon})\cap
W_{i})}{\rho(W_{i})}+\frac{\varepsilon}{2}\stackrel{{\scriptstyle\eqref{max_rho}}}{{<}}\varepsilon.$
(40)

Step 3. Let $x\in\mathbb{X}$ be an arbitrary initial condition. Then, there
exists $n_{2}=n_{2}(x)\in\mathbb{N}$ such that for all $n\geq n_{2}$ we have
${\mathbb{P}(\varphi_{q}^{n}(x)\not\in C_{\varepsilon})<\varepsilon}$.

Since $\varphi^{1}_{\mathcal{Q}}(x)$ is assumed to be finite, also
$\varphi^{r}_{\mathcal{Q}}(x):=\\{\varphi^{r}_{q}(x):q\in\mathcal{Q}\\}$ is
finite for any $r\in\mathbb{N}_{0}$. Let $n=k\ell+r$ for some
$k\in\mathbb{N}_{0}$ and $r\in\\{0,1,\ldots,\ell-1\\}$. Note that

$\mathbb{P}\left(\varphi_{q}^{\ell k+r}(x)\not\in
C_{\varepsilon}\right)=\sum_{y\in\varphi^{r}_{\mathcal{Q}}(x)}\mathbb{P}\left(\varphi^{\ell
k}_{q}(y)\not\in
C_{\varepsilon}\right)\cdot\mathbb{P}\left(\varphi^{r}_{q}(x)=y\right).$

Let $n_{1}(y)$ be given from Step 2 above, and choose
$k\geq\max\\{n_{1}(y):y\in\varphi^{r}_{\mathcal{Q}}(x)\\}$. Then

$\mathbb{P}(\varphi_{q}^{\ell k+r}(x)\not\in
C_{\varepsilon})\stackrel{{\scriptstyle\eqref{P_C_eps}}}{{<}}\varepsilon\cdot\sum_{y\in\varphi^{r}_{\mathcal{Q}}(x)}\mathbb{P}(\varphi^{r}_{q}(x)=y)\
=\varepsilon.$

So, the claim follows by taking $n_{2}(x):=\max\\{\ell\cdot
n_{1}(y):y\in\varphi^{r}_{\mathcal{Q}}(x),r\in\\{0,1,\ldots,\ell-1\\}\\}$.

Step 4. For each $K\subset\mathbb{N}_{0}$ compact (i.e. finite) there exists
$n_{0}=n_{0}(K)\in\mathbb{N}$ such that
$\mathbb{P}(\varphi_{q}^{n}(x)>z)<\varepsilon$ for all $n\geq n_{0}$ and $x\in
K$.

Notice that Step 3 assures that the claim holds when $K=\\{x\\}$ for any
$x\in\mathbb{N}_{0}$. If we consider any finite set
$K=\\{x_{1},x_{2},\ldots,x_{l}\\}$, then the result follows by taking
$n_{0}(K)=\max_{i=1,\ldots,l}n_{2}(x_{i})$. ∎

Note that, by virtue of Theorem 10, the existence of the weak attractor can be
derived purely by Markov chain arguments and is expected to occur in a large
class of RDS derived from reaction jump processes via the embedded Markov
chain approach. For instance, since chemical reaction networks are defined via
a finite set of reactions, $\varphi^{1}_{\mathcal{Q}}(x)$ is always finite for
any $x\in\mathbb{X}$. In particular, using Theorem 10 we can directly show
that the RDS corresponding to the embedded Markov chains of the birth-death
process and the Schlögl model admit a weak attractor, as presented in the next
sections.

###  4.3 Random periodic orbit of the birth-death chain

The concept of a random attractor makes sense only when considering an
invertible dynamical system on $\mathcal{Q_{+}}$. However, in the RDS
formulation of the embedded Markov chain (see Sec. 2.2), the shift map is not
invertible since for $q_{0}\neq\tilde{q}_{0}$ we have that

$\theta(q_{0},q_{1},q_{2},\ldots)=\theta(\tilde{q}_{0},q_{1},q_{2})=(q_{1},q_{2},\ldots).$

We can come around this inconvenience by redefining our noise space as

$\mathcal{Q}:=\\{q=(q_{n})_{n\in\mathbb{Z}}:q_{n}\in[0,1]\\},$

endowed with its Borel $\sigma$-algebra $\sigma(\mathcal{Q})$ and the bi-
infinite product measure $\lambda^{\mathbb{Z}}$. We redefine the shift map as

$\theta q=\theta(q_{n})_{n\in\mathbb{Z}}=(q_{n+1})_{n\in\mathbb{Z}}$

such that $(\theta^{-n}q)_{i}=q_{i-n}$, while the cocycle map $\varphi$
remains the same, see (11). Note that the synchronization results in Section 3
transfer immediately to the invertible setting since they have been formulated
independently from the past. By abuse of notation, we will from now on use
$\mathbb{P}=\lambda^{\mathbb{Z}}$.

####  4.3.1 Weak attraction

In this subsection we provide a full characterization of the weak random
attractor for the birth-death process given in Example 1.

######  Proposition 11.

The embedded Markov chain of the birth-death process admits a unique
stationary distribution. Therefore, the associated RDS $(\theta,\varphi)$
admits a (unique) weak attractor.

###### Proof.

Recall that the transition probabilities for the embedded Markov chain for the
birth-death process are given in (21).

From [8, p. 304], we know that any of the chains with transition probabilities
given by (21) or (32) admits a unique stationary distribution if and only if

$\zeta:=\sum_{x=1}^{\infty}\prod_{j=0}^{x-1}\frac{P_{1}(j)}{P_{-1}(j+1)}<\infty.$
(41)

We show that in both cases (41) holds. We consider the transition
probabilities (21), and for simplicity let
$\alpha:=\frac{\gamma_{2}}{\gamma_{1}}>0$. Then,

$\displaystyle\zeta$
$\displaystyle=\sum_{x=1}^{\infty}\prod_{j=0}^{x-1}\frac{1+\alpha(j+1)}{(1+\alpha
j)\alpha(j+1)}=\sum_{x=1}^{\infty}\prod_{j=0}^{x-1}\left(1+\frac{\alpha}{1+\alpha
j}\right)\cdot\frac{1}{\alpha(j+1)}$
$\displaystyle\quad\leq\sum_{x=1}^{\infty}\prod_{j=0}^{x-1}\left(\frac{1+\alpha}{\alpha}\right)\cdot\frac{1}{j+1}=\sum_{x=1}^{\infty}\frac{\beta^{x}}{x!},$
(42)

where $\beta=(1+\alpha)/\alpha$. Recall that by Stirling’s approximation we
have that $x!/(\sqrt{2\pi x}(\frac{x}{e})^{x})\rightarrow 1$ as
$x\rightarrow\infty$. Then, there exists $N\in\mathbb{N}$ such that for all
$x\geq N$ we have

$\frac{1}{x!}<\frac{2}{\sqrt{2\pi x}\left(\frac{x}{e}\right)^{x}}.$ (43)

From (4.3.1) and (43), and assuming without loss of generality that $N\geq
2\pi$, it follows that

$\zeta<S+\sqrt{2/\pi}\sum_{x=N}^{\infty}\beta^{x}\cdot\frac{1}{\sqrt{x}\left(\frac{x}{e}\right)^{x}}\leq
S+\sqrt{2/\pi}\sum_{x=N}^{\infty}\left(\frac{e\beta}{x}\right)^{x}$

for $S:=\sum_{x=1}^{N-1}\beta^{x}\cdot\frac{1}{x!}$. Hence, it suffices to
show that

$\sum_{x=N}^{\infty}\left(\frac{e\beta}{x}\right)^{x}<\infty.$

Assume again without loss of generality that $N$ is big enough so that
$\frac{e\beta}{x}<\frac{1}{e}$ for all $x\geq N$. Thus,

$\sum_{x=N}^{\infty}{\left(\frac{e\beta}{x}\right)^{x}}<\sum_{x=N}^{\infty}{e^{-x}}<\infty,$

from which we finally deduce $\zeta<\infty$. ∎

Figure 9 shows the stationary distribution of the one- and two-point motion
for the birth-death chain. The state space of the two-point motion splits up
into a transient area $\mathbb{X}^{2}\setminus\mathbb{D}$ and two
communication classes $\Delta$ and $\mathbb{D}\setminus\Delta=I_{-1}\cup
I_{1}$ [4] where $\Delta$ is the diagonal given in (17) and $\mathbb{D}$ is
the thick diagonal defined in (19). The respective stationary distributions
$\pi_{\Delta}$ and $\pi_{\mathbb{D}\setminus\Delta}$ are plotted in Figure
9(b)-(c).

\begin{overpic}[width=433.62pt]{Figures/BDP_RJP_stat_dist_bare.png}
\put(0.0,25.0){\scriptsize\begin{turn}{90.0} $\rho(x)$ \end{turn} }
\put(50.0,0.0){\scriptsize$x$} \end{overpic} (a)

\begin{overpic}[width=433.62pt]{Figures/BDP_RJP_2dim_stat_dist_even_bare.png}
\put(3.0,40.0){\scriptsize\begin{turn}{90.0} $y$ \end{turn} }
\put(43.0,2.0){\scriptsize$x$} \put(90.0,33.0){\scriptsize\begin{turn}{90.0}
$\pi_{\Delta}(x,y)$ \end{turn} } \end{overpic} (b)

\begin{overpic}[width=433.62pt]{Figures/BDP_RJP_2dim_stat_dist_odd_bare.png}
\put(3.0,40.0){\scriptsize\begin{turn}{90.0} $y$ \end{turn} }
\put(43.0,2.0){\scriptsize$x$} \put(93.0,33.0){\scriptsize\begin{turn}{90.0}
$\pi_{\mathbb{D}\setminus\Delta}(x,y)$ \end{turn} } \end{overpic} (c)

Figure 9: Stationary distributions of one- and two-point motion of the birth-
death chain. (a) Stationary distribution $\rho$ of the one-point motion and
(b), (c) stationary distributions $\pi_{\Delta}$ and
$\pi_{\mathbb{D}\setminus\Delta}$ of the two-point motion, respectively. All
stationary distributions are computed as the solution of their eigenvector
equation $\rho^{T}P=\rho^{T}$ for $P$ the corresponding transition matrix.
$\gamma_{1}=10,\gamma_{2}=1$.

In the following, we investigate the structure of the weak attractor by
formulating Proposition 12 and Corollary 13 as preparatory work for Section
4.3.2, where we will show that the weak attractor is in fact a strong pullback
attractor consisting of two random points that have distance one. The crucial
insight is the translation of the partial synchronization results from
Proposition 2 and Corollary 3 into the random periodic structure of the random
attractor. Depending on the type of partition in such a partial
synchronization, the following analysis may well be understood as a blueprint
for various forms of augmented Markov chains with more complicated structure.

######  Proposition 12.

The weak attractor $q\mapsto A_{q}$ of the birth-death chain has two points on
each fiber $\mathbb{P}$-a.s., that is

$\mathbb{P}(\\#A_{q}=2)=1.$

Moreover, $A_{q}\cap W_{i}\neq\emptyset$ for $i=0,1$ $\mathbb{P}$-a.s., where
$W_{i}$ are given in Corollary 3.

###### Proof.

We split the proof in two parts. We first show that $A_{q}$ has at least two
points $\mathbb{P}$-a.s., and then we conclude that $\mathbb{P}(\\#A_{q}\geq
3)=0$.

Step 1. $A_{q}\cap W_{i}\neq\emptyset$ for $i=0,1$ $\mathbb{P}$-a.s.

For each $x\in\mathbb{N}_{0}$ consider the pullback limit

$a_{x}^{+}(q):=\limsup_{n\rightarrow\infty}\varphi^{2n}_{\theta^{-2n}q}(x).$

Let $(n_{j})_{j=0}^{\infty}$ be a strictly increasing sequence of natural
numbers such that $\varphi^{2n_{j}}_{\theta^{-2n_{j}}q}(x)\rightarrow
a_{x}^{+}(q)$ for $j\to\infty$. Since $A_{q}$ is a weak attractor we have

$\lim_{j\rightarrow\infty}\mathbb{P}\left(d\left(\varphi^{2n_{j}}_{\theta^{-2n_{j}}q}(x),A_{q}\right)\geq
1\right)=0.$

From the last limit we consider a subsequence $n_{j_{k}}$ such that
$\varphi^{2n_{j_{k}}}_{\theta^{-2n_{j_{k}}}q}(x)\in A_{q}$ for all
$k\in\mathbb{N}_{0}$ $\mathbb{P}$-a.s. We conclude that $a_{x}^{+}(q)\in
A_{q}$ for $\mathbb{P}$-a.e. $q$, which in particular implies that
$a_{x}^{+}<\infty$ $\mathbb{P}$-a.s. Furthermore, since $a_{x}^{+}(q)\in
W_{i}$ if and only if $x\in W_{i}$, the claim follows by taking $x=0,1$, for
instance.

Step 2. $\mathbb{P}(\\#A_{q}\geq 3)=0$.

Since $\theta$ is $\mathbb{P}$-invariant, we have that for any
$n\in\mathbb{N}$

$\mathbb{P}(\\#A_{q}\geq 3)=\mathbb{P}(\\#A_{\theta^{n}q}\geq
3)=\mathbb{P}(\\#\varphi^{n}_{q}(A_{q})\geq 3),$

where the last equality follows by the invariance of $A_{q}$. We now partition

$\varphi^{n}_{q}(A_{q})=\varphi^{n}_{q}(A_{q}\cap
W_{0})\cup\varphi^{n}_{q}(A_{q}\cap W_{1}).$

Since it follows from Step 1 that for all $n\in\mathbb{N}$ and $i\in\\{0,1\\}$
the sets $\varphi^{n}_{q}(A_{q}\cap W_{i})$ contain at least one element, we
can combine these observations to give the bound

$\mathbb{P}(\\#A_{q}\geq 3)\leq\mathbb{P}(\\#\varphi^{n}_{q}(A_{q}\cap
W_{0})\geq 2)+\mathbb{P}(\\#\varphi^{n}_{q}(A_{q}\cap W_{1})\geq 2).$

By taking the limit $n\rightarrow\infty$, Proposition 6 implies that the
right-hand side tends to $0$ by taking $K^{i}_{q}=A_{q}\cap W_{i}$. The result
follows. ∎

######  Corollary 13.

For each $x\in\mathbb{N}_{0}$

$a_{x}(q):=\lim_{n\rightarrow\infty}\varphi^{2n}_{\theta^{-2n}q}(x)$ (44)

exists $\mathbb{P}$-a.s.. Furthermore, for $x,y\in W_{i}$, $i=0,1$, we have
that $a_{x}=a_{y}$ $\mathbb{P}$-a.s.. Conversely, if $x\in W_{0}$ and $y\in
W_{1}$, then $a_{x}\neq a_{y}$ $\mathbb{P}$-a.s..

###### Proof.

For each $x\in\mathbb{N}_{0}$ let

$a_{x}^{-}(q):=\liminf_{n\rightarrow\infty}\varphi^{2n}_{\theta^{-2n}q}(x).$

Analogously to Step 1 in the proof of Proposition 12, we obtain
$a_{x}^{-}(q)\in A_{q}$ for almost all $q\in\mathcal{Q}$. Since it has the
same parity as $a_{x}^{+}(q)$ and $\\#A_{q}\cap W_{i}=1$, $i=0,1$, we conclude
that $a_{x}^{+}=a_{x}^{-}$ holds with full probability and, hence, the limit
$a_{x}\in A_{q}$ exists almost surely.

Using again that $\\#A_{q}\cap W_{i}=1$, $i=0,1$, we derive that $a_{x}$ and
$a_{y}$ are identical (or different, respectively) in a full measure set when
$x$ and $y$ are of the same parity (or of different parities, respectively). ∎

####  4.3.2 Strong pullback attraction to a random period orbit

We can now give a full characterization of the pullback attractor in the next
theorem. We already know from Proposition 12 and Corollary 13 that the weak
attractor $A_{q}$ consists almost surely of the two distinct random points
$a_{0}(q)$ and $a_{1}(q)$, as given in (44). Now, we identify the strong
pullback structure of this attractor.

######  Theorem 14.

The weak attractor $A_{q}=\\{a_{0}(q),a_{1}(q)\\}$ for the birth-death chain

  * (i)

is a random periodic orbit of period 2,

  * (ii)

is a pullback and a forward attractor, and

  * (iii)

satisfies $\left|a_{0}(q)-a_{1}(q)\right|=1$ for $\mathbb{P}$-almost all
$q\in\mathcal{Q}$.

###### Proof.

At first, we use the cocycle property in order to show item (i), that is

$\varphi_{q}^{1}(a_{0}(q))=a_{1}(\theta
q)\qquad\text{and}\qquad\varphi_{q}^{1}(a_{1}(q))=a_{0}(\theta q)$ (45)

is satisfied $\mathbb{P}$-a.s. Indeed, let $i=0,1$, $q\in\mathcal{Q}$ be
fixed. Then,

$\varphi_{q}^{1}(a_{i}(q))=\lim_{n\rightarrow\infty}\varphi^{2n+1}_{\theta^{-2n}q}(i)=\lim_{n\rightarrow\infty}\varphi^{2n}_{\theta^{-2n}\circ\theta
q}\left(\varphi_{q}^{1}(i)\right)=a_{i+1(\bmod 2)}(\theta q),$

where the last equality follows directly from Corollary 13 and the fact that
$\varphi_{q}^{1}(i)\in W_{i+1(\bmod 2)}$.

In order to show (ii), let $q\in\mathcal{Q}$. It follows from the definition
(44) that, for any $x\in W_{i}$ with $i=0,1$, there exists
$N_{0}=N_{0}(x)\in\mathbb{N}_{0}$ such that for all $n\geq N_{0}$ we have
$\varphi_{\theta^{-2n}q}^{2n}(x)=a_{i}(q)\in A(q)$. On the other hand, by the
cocycle property and (45) it follows that

$\displaystyle\lim_{n\rightarrow\infty}\varphi_{\theta^{-2n-1}}^{2n+1}(x)$
$\displaystyle=\lim_{n\rightarrow\infty}\varphi^{1}_{\theta^{-1}q}\left(\varphi^{2n}_{\theta^{-2n}\theta^{-1}q}(x)\right)$
$\displaystyle=\varphi^{1}_{\theta^{-1}q}(a_{i}(\theta^{-1}q))=a_{i+1(\bmod
2)}(q).$

Hence, there is $N_{1}=N_{1}(x)\in\mathbb{N}$ such that for all $n\geq N_{1}$
we have that $\varphi_{\theta^{-2n-1}q}^{2n+1}(x)\in A(q)$. Combining both
parts, we obtain $\varphi_{\theta^{-n}q}^{n}(x)\in A(q)$ for all $n\geq
N(x):=\max\\{2N_{0}(x),2N_{1}(x)+1\\}$. Recall that point pullback attractors
are (set) pullback attractors due to the state space being discrete, and thus
$A$ is a pullback attractor. $A$ is also a forward attractor due to Theorem 8.

Last, we prove (iii). Consider the sets

$\Omega_{n}:=\left\\{q\in\mathcal{Q}:\varphi_{\theta^{-2k}q}^{2k}(0)=a_{0}(q),\
\varphi_{\theta^{-2k}q}^{2k}(1)=a_{1}(q)\ \forall k\geq n\right\\}.$

Since $\Omega_{n}\subset\Omega_{n+1}$ for all $n\in\mathbb{N}$, and since
$\mathbb{P}(\bigcup_{n\in\mathbb{N}}\Omega_{n})=1$, for each $\varepsilon>0$
there is $N$ such that for all $k\geq N$ sufficiently large

$\displaystyle\mathbb{P}\left(\left|a_{0}(q)-a_{1}(q)\right|\geq 2\right)$
$\displaystyle\leq\mathbb{P}\left(\\{q\in\mathcal{Q}:\left|a_{0}(q)-a_{1}(q)\right|\geq
2\\}\ \cap\Omega_{N}\right)+\frac{\varepsilon}{2}$
$\displaystyle=\mathbb{P}\left(\left|\varphi_{\theta^{-2k}q}^{2k}(0)-\varphi_{\theta^{-2k}q}^{2k}(1)\right|\geq
2\right)+\frac{\varepsilon}{2}$
$\displaystyle=\mathbb{P}\left(\left|\varphi_{q}^{2k}(0)-\varphi_{q}^{2k}(1)\right|\geq
2\right)+\frac{\varepsilon}{2}<\varepsilon,$

due to Proposition 2. Since $\varepsilon$ was arbitrarily small, the result
follows. ∎

####  4.3.3 Sample measures supported on the attractor

In this subsection we briefly describe the statistical importance of the
attractor $A$ in terms of the invariant measure for the skew-product map
$\Theta:\mathcal{Q}\times\mathbb{X}\rightarrow\mathcal{Q}\times\mathbb{X}$
given by

$\Theta(q,x):=(\theta q,\varphi_{q}(x)).$

Denoting by $T^{*}\mu$ the push forward of a measure $\mu$ by a map $T$, i.e.
$T^{*}\mu(\cdot)=\mu(T^{-1}(\cdot))$, we adopt the classical definition of an
invariant measure for the RDS (see e.g. [3, Definition 1.4.1]): A probability
measure $\mu$ on $\mathcal{Q}\times\mathbb{X}$ is invariant for the random
dynamical system $(\theta,\varphi)$ if

  1. (i)

$\Theta_{t}^{*}\mu=\mu$ for all $t\in\mathbb{\mathbb{missing}}{N}_{0}$ ,

  2. (ii)

the marginal of $\mu$ on $\mathcal{Q}$ is $\mathbb{P}$, i.e. $\mu$ can be
factorized uniquely into

$\mu(\mathrm{d}q,x)=\mu_{q}(x)\mathbb{P}(\mathrm{d}q),$

where $q\mapsto\mu_{q}$ is the sample measure (or disintegration) on
$\mathbb{X}$, i.e., $\mu_{q}$ is almost surely a probability measure on
$\mathbb{X}$ and $q\mapsto\mu_{q}(B)$ is measurable for all
$B\subset\mathbb{X}$.

In particular note that, since $\mathbb{P}$ is given, the sample measures
$\mu_{q}$ completely determine such an invariant measure. A specific form of
such invariant measures are _Markov measures_ , characterized by the sample
measures being measurable with respect to the _past_ : in our setting, this
means that the $\mu_{q}$ only depend on $q_{n},n<0$ (cf. e.g. [24] or [20]).

The theory of random dynamical systems now gives us the following result on
the unique invariant measure for the RDS at hand, relating it to the
stationary distribution of the Markov chain:

######  Proposition 15.

The RDS of the birth-death chain possesses a unique invariant Markov measure
with sample measures

$\mu_{q}=\frac{1}{2}\delta_{a_{0}(q)}+\frac{1}{2}\delta_{a_{1}(q)},$

such that $\mathbb{E}[\mu_{q}]=\rho$, where $\rho$ is the unique stationary
distribution from Proposition 11.

###### Proof.

By [6, Proposition 4.5], we know that there exists a Markov measure $\mu$ such
that $\mu_{q}(A(q))=1$ almost surely, where $A(q)=\\{a_{0}(q),a_{1}(q)\\}$ is
the attractor from Theorem 14. Its uniqueness and the fact that
$\mathbb{E}[\mu_{q}]=\rho$ follow from the celebrated correspondence theorem,
also called _Ledrappier-LeJan-Crauel Theorem_ (see [24, Proposition 1.2.3] for
a version that suffices for our situation and [20, Theorem 4.2.9] for the more
general situation).

Using [23] (see also [11, Lemma 2.19]) we can directly infer that either
$\mu_{q}=\frac{1}{2}\delta_{a_{0}(q)}+\frac{1}{2}\delta_{a_{1}(q)}$ almost
surely or $\mu_{q}=\delta_{a_{i}(q)}$ almost surely for $i=0$ or $i=1$ fixed.
The latter case can now be excluded by combining (45) and the invariance
property $(\varphi_{q}^{n})^{*}\mu_{q}=\mu_{\theta^{n}q}$ [20, Proposition
1.3.27]. ∎

###  4.4 Attractor structure for the Schlögl model

For the extended reaction network of Example 2, we show the existence of a
weak attractor for the corresponding embedded Markov chain via its possession
of a unique stationary distribution. This is formulated in the following
proposition. As for the structure of the weak attractor, we provide a
conjecture based on numerical experiments whose proof is left as an open
problem for future work.

######  Proposition 16.

The embedded Markov chain of the Schlögl model admits a unique stationary
distribution. Therefore, the associated RDS admits a (unique) weak attractor.

###### Proof.

Analogously to Proposition 11, we show that the quantity $\zeta$ in (41) is
finite. Indeed,

$\prod_{j=0}^{x-1}\frac{P_{1}(j)}{P_{-1}(j+1)}=\prod_{j=0}^{x-1}\left(\frac{\gamma_{1}+\gamma_{3}j(j-1)}{\gamma_{2}(j+1)+\gamma_{4}j(j+1)(j-1)}\right)\cdot\left(\frac{\mu(j+1)}{\mu(j)}\right)=$
$=\frac{1}{x!}\prod_{j=0}^{x-1}\left(\frac{\gamma_{1}+\gamma_{3}j(j-1)}{\gamma_{2}+\gamma_{4}j(j-1)}\right)\cdot\left(\frac{\mu(j+1)}{\mu(j)}\right).$

Both terms inside the product are bounded for all $j\in\mathbb{N}_{0}$. Thus,
there exists $C>0$ such that

$\zeta\leq\sum_{x=1}^{\infty}\frac{1}{x!}C^{x}<\infty.$

The finiteness of the right-hand side follows from Stirling’s approximation
similarly as in Proposition 11. ∎

The stationary distribution of the Schlögl model is depicted in Figure 10.

Numerical explorations suggest the following structure of the weak attractor:
similarly to the simple birth-death chain, there is partial synchronization
with respect to a partition of $\mathbb{N}_{0}$ into odd and even numbers, as
exemplified in Figure 7(a). This suggests that the weak attractor consists
again of two random points that most likely form a random periodic orbit.
However, as we see in Figure 7(b), the distance of these two random points
will not be one but probably larger; maybe even depending on the random
realization. We emphasize that this difference to the simple birth-death chain
in terms of the non-synchronizing trajectories reflects the distinction
between monostability (the thick diagonal in the two-point motion is
absorbing) and bistability (the thick diagonal in the two-point motion has
repelling parts) as seen in the large-volume limiting ODEs (cf. Figure 1).
This fact is also mirrored by the respective stationary distributions of the
one- and two-point motion, as illustrated in Figures 9 and 10.

A proof of the associated structure of the random attractor is much more
involved than in the simple birth-death case and will be left for future work.

\begin{overpic}[width=433.62pt]{Figures/BDPS_RJP_stat_dist_alt2_bare.png}
\put(-2.0,25.0){\scriptsize\begin{turn}{90.0} $\rho(x)$ \end{turn} }
\put(50.0,0.0){\scriptsize$x$} \end{overpic} (a)

\begin{overpic}[width=433.62pt]{Figures/BDPs_RJP_2dim_stat_dist_even_bare.png}
\put(3.0,40.0){\scriptsize\begin{turn}{90.0} $y$ \end{turn} }
\put(43.0,2.0){\scriptsize$x$} \put(93.0,33.0){\scriptsize\begin{turn}{90.0}
$\pi_{\Delta}(x,y)$ \end{turn} } \end{overpic} (b)

\begin{overpic}[width=433.62pt]{Figures/BDPS_RJP_2dim_stat_dist_odd_bare.png}
\put(3.0,40.0){\scriptsize\begin{turn}{90.0} $y$ \end{turn} }
\put(43.0,2.0){\scriptsize$x$} \put(95.0,33.0){\scriptsize\begin{turn}{90.0}
$\pi_{S}(x,y)$ \end{turn} } \end{overpic} (c)

\begin{overpic}[width=433.62pt]{Figures/BDPS_RJP_2dim_stat_dist_odd2_bare.png}
\put(3.0,40.0){\scriptsize\begin{turn}{90.0} $y$ \end{turn} }
\put(45.0,2.0){\scriptsize$x$} \put(93.0,33.0){\scriptsize\begin{turn}{90.0}
$\pi_{S}(x,y)$ \end{turn} } \end{overpic} (d)

Figure 10: Stationary distributions of one- and two-point motion of the
Schlögl model (embedded Markov chain). (a) Stationary distribution $\rho$ of
the one-point motion and (b), (c) stationary distributions $\pi_{\Delta}$ and
$\pi_{S}$ of the two-point motion, respectively. Here, the index
$S\subset\mathbb{X}^{2}\setminus\Delta$ refers to the support of the
stationary distribution for the two-point motion starting with odd distance
between the two coordinates $x_{0},y_{0}$ of the initial state. (d) the same
as (c) only on a logarithmic scale. $\gamma_{1}=6$,
$\gamma_{2}=3.5,\gamma_{3}=0.4,\gamma_{4}=0.0105$.

##  5 Conclusion

We have introduced the phenomenon of (partial) time-shifted synchronization
for reaction jump processes, using their description via the augmented and
embedded Markov chain whose properties as a random dynamical system can be
specified through the structure of the corresponding random attractor. As a
first example we have given a full proof of partial synchronization for the
birth-death process, finding the random attractor of the embedded chain to be
a random periodic orbit of period $2$. We have demonstrated that for
extensions of this basic example, such as the Schlögl model, one may expect a
similar structure; however, there is an apparent difference in the relation of
the two random points.

Depending on the dimensionality of the system, i.e. the number of different
species, or other types of reaction rates, e.g. Michaelis-Menten, we expect
various forms of (partial) synchronization and random periodic behavior for
general reaction systems and leave it as a research direction for the future
to work towards a categorization of such processes and their corresponding
chains in terms of random dynamical systems theory. The main goal of this
paper has been to relate synchronization phenomena for the different
formalizations of the chemical reaction process and to give a first complete
and rigorous description of a classical example. The strategy of understanding
the two-point motion, establishing the attractor via a unique stationary
distribution with sufficient decay and then specifying the structure of the
attractor and the statistical equilibrium supported there, based on the two-
point motion analysis, may well be generalizable. In particular, one may
consider time-shifted synchronization for general Markov jump processes, not
necessarily given by reaction networks, via the RDS description of the related
space-time Markov chains. Moreover, an intriguing point of more detailed
investigation will concern the quantification and statistics of the delay
times found for time-shifted synchronization, which may be of high interest
also for the applied side of chemical reaction processes.

Furthermore, our work has brought up additional questions that remain open, to
our knowledge. Can one find an example coming from a chemical reaction network
with no (partial) synchronization at all, i.e. where each synchronization
class is a singleton? May one describe bifurcations of the attractor, for
example in an easy model such as Schlögl’s, via variation of the parameters?
How are attractors of the described processes related to the attractors of the
corresponding volume-scaled systems, i.e. the Langevin SDEs or the reaction
rate ODEs? Additionally, from the RDS point of view it will also be intriguing
to give general criteria for weak attractors being (strong) pullback
attractors in discrete state spaces. In summary, we see this work as a first
step towards a deeper structural understanding of reaction jump processes via
RDS theory and, conversely, a motivation for a broader understanding of random
attractors withing the dichotomy between the discrete and the continuous.

### Acknowledgements

We acknowledge the support of Deutsche Forschungsgemeinschaft (DFG) through
CRC 1114 and under Germany’s Excellence Strategy – The Berlin Mathematics
Research Center MATH+ (EXC-2046/1, project ID 390685689). M. E. additionally
thanks the DFG-funded SPP 2298 for supporting his research. G. O.-M. also
thanks FU Berlin for a 3-month Forschungsstipendium. The authors gratefully
acknowledge Dennis Chemnitz for fruitful discussions.

## References

  * Anderson and Kurtz [2011] D. F. Anderson and T. G. Kurtz.  Continuous time Markov chain models for chemical reaction networks.  In H. Koeppl, G. Setti, M. di Bernardo, and D. Densmore, editors, _Design and analysis of biomolecular circuits_ , pages 3–42. Springer, New York, NY, 2011.  doi: 10.1007/978-1-4419-6766-4˙1.
  * Anderson and Kurtz [2015] D. F. Anderson and T. G. Kurtz.  _Stochastic analysis of biochemical systems_ , volume 674.  Springer, 2015.  doi: 10.1007/978-3-319-16895-1.
  * Arnold [1998] L. Arnold.  _Random Dynamical Systems_.  Springer, Berlin, 1998.  doi: 10.1007/BFb0095238.
  * Baxendale [1991] P. H. Baxendale.  Statistical equilibrium and two-point motion for a stochastic flow of diffeomorphisms.  In K. S. Alexander and J. C. Watkins, editors, _Spatial stochastic processes_ , pages 189–218. Springer, 1991.  doi: 10.1007/978-1-4612-0451-0˙9.
  * Crauel [2002] H. Crauel.  _Random probability measures on Polish spaces, vol. 11 of Stochastics Monographs_.  Taylor & Francis, 2002.
  * Crauel and Flandoli [1994] H. Crauel and F. Flandoli.  Attractors for random dynamical systems.  _Probab. Theory Relat. Fields_ , 100(3):365–393, 1994.  doi: 10.1007/BF01193705.
  * Crauel and Kloeden [2015] H. Crauel and P. Kloeden.  Nonautonomous and random attractors.  _Jahresbericht der Deutschen Mathematiker-Vereinigung_ , 117:173–206, 06 2015.  doi: 10.1365/s13291-015-0115-0.
  * Durrett [2010] R. Durrett.  _Probability: Theory and Examples_ , volume 49 of _Cambridge Series in Statistical and Probabilistic Mathematics_.  Cambridge University Press, 4th edition, 2010.
  * Endres [2017] R. Endres.  Entropy production selects nonequilibrium states in multistable systems.  _Sci. Rep._ , 7:14437, 10 2017.  doi: 10.1038/s41598-017-14485-8.
  * Engel and Kuehn [2021] M. Engel and C. Kuehn.  A random dynamical systems perspective on isochronicity for stochastic oscillations.  _Comm. Math. Phys._ , 386(3):1603–1641, 2021\.  ISSN 0010-3616.  doi: 10.1007/s00220-021-04077-z.
  * Flandoli et al. [2017] F. Flandoli, B. Gess, and M. Scheutzow.  Synchronization by noise.  _Probab. Theory Relat. Fields_ , 168(3):511–556, 2017.  doi: 10.1007/s00440-016-0716-2.
  * Gillespie [1976] D. T. Gillespie.  A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.  _J. Comput. Phys._ , 22(4):403–434, 1976.  doi: 10.1016/0021-9991(76)90041-3.
  * Gillespie [1977] D. T. Gillespie.  Exact stochastic simulation of coupled chemical reactions.  _J. Phys. Chem._ , 81(25):2340–2361, 1977.  doi: 10.1021/j100540a008.
  * Gillespie [1992] D. T. Gillespie.  A rigorous derivation of the chemical master equation.  _Physica A_ , 188(1-3):404–425, 1992.  doi: 10.1016/0378-4371(92)90283-V.
  * Gillespie [2000] D. T. Gillespie.  The chemical Langevin equation.  _J. Chem. Phys._ , 113(1):297–306, 2000.  doi: 10.1063/1.481811.
  * Huang et al. [2020] W. Huang, H. Qian, S. Wang, F. X.-F. Ye, and Y. Yi.  Synchronization in discrete-time, discrete-state random dynamical systems.  _SIAM J. Appl. Dyn. Syst._ , 19(1):233–251, 01 2020.  doi: 10.1137/19M1244883.
  * Jahnke [2011] T. Jahnke.  On reduced models for the chemical master equation.  _Multiscale Modeling & Simulation_, 9(4):1646–1676, 2011.  doi: 10.1137/110821500.
  * Kifer [1986] Y. Kifer.  _Ergodic Theory of Random Transformations_.  Birkhäuser Boston, 1 1986.  ISBN 978-1-4684-9177-7.  doi: 10.2307/2288883.
  * Ko [1991] M. S. Ko.  A stochastic model for gene induction.  _J. Theor. Biol._ , 153(2):181–194, 1991.  doi: 10.1016/S0022-5193(05)80421-7.
  * Kuksin and Shirikyan [2012] S. Kuksin and A. Shirikyan.  _Mathematics of two-dimensional turbulence_ , volume 194 of _Cambridge Tracts in Mathematics_.  Cambridge University Press, 2012.  ISBN 978-1-107-02282-9.  doi: 10.1017/CBO9781139137119.
  * Kurtz [1970] T. G. Kurtz.  Solutions of ordinary differential equations as limits of pure jump Markov processes.  _J. Appl. Probab._ , 7(1):49–58, 1970.  doi: 10.2307/3212147.
  * Kurtz [1972] T. G. Kurtz.  The relationship between stochastic and deterministic models for chemical reactions.  _J. Chem. Phys._ , 57(7):2976–2978, 1972.  doi: 10.1063/1.1678692.
  * Le Jan [1987] Y. Le Jan.  Équilibre statistique pour les produits de difféomorphismes aléatoires indépendants.  _Ann. Inst. H. Poincaré Probab. Statist._ , 23(1):111–120, 1987.
  * Ledrappier and Young [1994] F. Ledrappier and L.-S. Young.  Entropy formula for random transformations.  _Probab. Theory Relat. Fields_ , 80(3):217–240, 1994.  doi: 10.1007/BF00356103.
  * Matheson et al. [1975] I. Matheson, D. F. Walls, and C. W. Gardiner.  Stochastic models of firstorder nonequilibrium phase transitions in chemical reactions.  _J Stat Phys_ , 12(1):21–34, 08 1975.  doi: 10.1007/BF01024182.
  * McQuarrie [1967] D. A. McQuarrie.  Stochastic approach to chemical kinetics.  _J. Appl. Probab._ , 4(3):413–478, 1967.  doi: 10.2307/3212214.
  * Menz et al. [2012] S. Menz, J. C. Latorre, C. Schutte, and W. Huisinga.  Hybrid stochastic–deterministic solution of the chemical master equation.  _Multiscale Modeling & Simulation_, 10(4):1232–1262, 2012.  doi: 10.1137/110825716.
  * Newman [2018] J. Newman.  Necessary and sufficient conditions for stable synchronization in random dynamical systems.  _Ergod. Theory Dyn. Syst._ , 38(5):1857–1875, 2018.  ISSN 0143-3857.  doi: 10.1017/etds.2016.109.
  * Newman [2020] J. Newman.  Synchronisation of almost all trajectories of a random dynamical system.  _Discrete Contin. Dyn. Syst._ , 40(7):4163–4177, 2020.  ISSN 1078-0947.  doi: 10.3934/dcds.2020176.
  * Schlögl [1972] F. Schlögl.  Chemical reaction models for non-equilibrium phase transitions.  _Z. Physik_ , 253(2):147–161, 1972.  doi: 10.1007/BF01379769.
  * Sikorski et al. [2021] A. Sikorski, M. Weber, and C. Schütte.  The augmented jump chain.  _Adv. Theory Simul._ , 4:2000274, 03 2021.  doi: 10.1002/adts.202000274.
  * Wilkinson [2019] D. J. Wilkinson.  _Stochastic modelling for systems biology_.  Chapman and Hall/CRC, 3rd edition, 2019.  doi: 10.1201/9781351000918.
  * Winkelmann and Schütte [2017] S. Winkelmann and C. Schütte.  Hybrid models for chemical reaction networks: Multiscale theory and application to gene regulatory systems.  _J. Chem. Phys._ , 147(11):114115, 2017.  doi: 10.1063/1.4986560.
  * Winkelmann and Schütte [2020] S. Winkelmann and C. Schütte.  _Stochastic Dynamics in Computational Biology_.  Springer, 2020.  doi: 10.1007/978-3-030-62387-6.
  * Ye and Qian [2019] F. X.-F. Ye and H. Qian.  Stochastic dynamics II: Finite random dynamical systems, linear representation, and entropy production.  _Discrete Contin. Dyn. Syst. Ser. B_ , 22(8):4341–4366, 04 2019.  doi: 10.3934/dcdsb.2019122.
  * Ye et al. [2016] F. X.-F. Ye, Y. Wang, and H. Qian.  Stochastic dynamics: Markov chains and random transformations.  _Discrete Contin. Dyn. Syst. Ser. B_ , 21(7):2337–2361, 08 2016.  doi: 10.3934/dcdsb.2016050.
  * Zeiser et al. [2008] S. Zeiser, U. Franz, O. Wittich, and V. Liebscher.  Simulation of genetic networks modelled by piecewise deterministic Markov processes.  _IET Syst Biol._ , 2(3):113–135, 2008.  doi: 10.1049/iet-syb:20070045.
  * Zhao and Zheng [2009] H. Zhao and Z.-H. Zheng.  Random periodic solutions of random dynamical systems.  _J. Differential Equations_ , 246(5):2020–2038, 2009.  doi: 10.1016/j.jde.2008.10.011.