# Noise in random Boolean networks

Tiago P. Peixoto<EMAIL_ADDRESS>Barbara Drossel<EMAIL_ADDRESS>darmstadt.de Institut für Festkörperphysik, TU Darmstadt, Hochschulstrasse 6,
64289 Darmstadt, Germany

###### Abstract

We investigate the effect of noise on Random Boolean Networks. Noise is
implemented as a probability $p$ that a node does not obey its deterministic
update rule. We define two order parameters, the long-time average of the
Hamming distance between a network with and without noise, and the average
frozenness, which is a measure of the extent to which a node prefers one of
the two Boolean states. We evaluate both order parameters as function of the
noise strength, finding a smooth transition from deterministic ($p=0$) to
fully stochastic ($p=1/2$) dynamics for networks with $K\leq 2$, and a first
order transition at $p=0$ for $K>2$. Most of the results obtained by computer
simulation are also derived analytically. The average Hamming distance can be
evaluated using the annealed approximation. In order to obtain the
distribution of frozenness as function of the noise strength, more
sophisticated self-consistent calculations had to be performed. This
distribution is a collection of delta peaks for $K=1$, and it has a fractal
sructure for $K>1$, approaching a continuous distribution in the limit $K\gg
1$.

###### pacs:

89.75.Da,05.65.+b,91.30.Dk,91.30.Px

##  I Introduction

Random Boolean Networks (RBNs) Kauffman (1969); Drossel (2008) have been used
as a simple model for a variety of dynamical systems consisting of interacting
units, such as neural networks Rosen-Zvi et al. (2001), social networks
Moreira et al. (2004) and, more prominently, gene regulatory networks
Lagomarsino et al. (2005); Kauffman (1969). RBNs are composed of Boolean nodes
that are coupled to each other. In the case of gene regulatory networks, the
Boolean state is a step-function approach to the expression level of a
particular gene. Despite this loss of detail, the most important features of
gene regulatory processes are still captured in many cases, since they should
not depend on biochemical details, but on the desired sequence of events in
the cell Bornholdt (2005).

So far, the dynamics of RBNs have mostly been studied using deterministic
update rules. The dynamics of such models are non-ergodic, with periodic
attractor trajectories in state space. Once the system has reached an
attractor, it remains there. Another important property of RBNs is a phase
transition, which occurs when the number $K$ of inputs per node is changed.
For unbiased networks, the dynamics exhibit a frozen phase at $K=1$, where
local perturbations die out quickly and most attractors are fixed points, and
a “chaotic” phase at $K>2$, where perturbations increase exponentially fast
and attractors have very long periods. At the boundary $K=2$ between those two
phases are the so-called “critical” networks, where perturbations increase
algebraically with time. Originally, it was suggested by Kauffmann Kauffman
(1969) that such critical networks are best suited to model real systems,
which are supposedly poised “at the edge of chaos”. In the meantime, there is
agreement that RBNs of all three types have only limited validity when applied
to real systems.

Real networks usually have some level of stochastic behaviour, and for this
reason several authors have investigated RBNs under the influence of
stochasticity. For instance, in Greil and Drossel (2005) the nodes of RBNs
were updated in a completely random order. This update method preserves the
non-ergodicity of the system, and it is still possible to identify distinct
attractors. Attractors are in this case defined as sets of states all of which
are visited for a non-vanishing proportion of time during the same trajectory.
The stochastic update sequence vastly reduces the number of attractors of
critical RBNs, which becomes a power law as function of the network size.
Similar results are obtained when the update sequence deviates only slightly
from a synchronous update Klemm and Bornholdt (2005). Such a power law was for
a long time falsely believed to occur in deterministic RBNs Kauffman (1969);
Samuelsson and Troein (2003); Drossel (2005); Drossel et al. (2005).

Instead of introducing stochasticity into the update times, other authors
introduce it into the update functions. In Shmulevich et al. (2002),
probabilistic Boolean functions are used, where a set of several Boolean
functions is assigned to each node, and at each time step one of these is
chosen randomly with a given probability. According to Shmulevich et al.
(2002), this model is more realistic than models with a purely deterministic
update scheme.

However, the most important way of introducing noise into a RBN is in form of
a “temperature”, leading often to ergodic behavior. The effect of thermal
noise on Ising spins on a network was studied in Indekeu (2004); Aleksiejuk et
al. (2002), where a “ferromagnetic” transition from the ordered to the
disordered phases was observed at a critical noise strength value. In the
language of gene regulatory networks, a temperature manifests itself as
fluctuations in the protein concentrations, so that a gene may not always be
turned on or off, given the same expression state of the other genes McAdams
and Arkin (1997). This effect can be included into models by allowing a
deviation from the deterministic update rule with a certain probability. In
Fretter and Drossel (2008), for instance, a subset of nodes were perturbed in
this way (this corresponds to turning on the temperature for a short time
interval), and the response of the dynamics to this perturbation was
evaluated, giving information about the basin structure of the system. Miranda
et al. Miranda and Parga (1989) studied the effect of a permanently acting
temperature by introducing a fixed probability $p$ that the state of a node
becomes the opposite of what it should be according to the deterministic
update rule. They evaluated the average crossing time between trajectories in
state space which started from different initial states as function of noise
strength Miranda and Parga (1989); Golinelli and Derrida (1989); Qu et al.
(2002). By sampling the entire state space of small networks ($N\leq 20$), it
was found that the “barriers”, which correspond to the attractor basin
boundaries, can be crossed with non-vanishing probability when $p>0$, although
the characteristic times may be large. This means that the system is always
ergodic. This type of noise has also been studied for Boolean networks with
threshold functions, corresponding to a majority update rule Huepe and Aldana-
González (2002). This system undergoes a second order phase transition at a
critical noise strength from an ordered dynamical phase, where all nodes
assume the same value for the majority of time, to a disordered phase where
nodes assume both states equally often.

In this work, we investigate the effect of ongoing stochastic noise on RBNs.
Following Miranda and Parga (1989); Golinelli and Derrida (1989); Qu et al.
(2002), noise strength is tuned via a probability $p$ that a node does not
obey its deterministic update rule. We monitor the transition from fully
deterministic dynamics ($p=0$) to purely stochastic dynamics ($p=1/2$) as the
noise strength is varied. Differently from Miranda and Parga (1989); Golinelli
and Derrida (1989); Qu et al. (2002), we are interested in the behaviour of
the networks in the limit of large system size, where it is impossible to
explore large parts of the state space. In order to characterize the
transition from zero to infinite temperature, we define two order parameters:
the long-time average of the Hamming distance between a network with and
without noise, and the average frozenness, which is a measure of the extent to
which a node prefers one of the two Boolean states. We find, both analytically
and numerically, that this transition for the Hamming distance is continuous
for $K\leq 2$, and discontinuous at $p=0$ for $K>2$, when the Hamming distance
is considered. This distinction is a direct consequence of the phase
transition from frozen to chaotic dynamics in the deterministic model. The
frozenness shows a smooth transition for all values of $K$. The distribution
of frozenness shows a surprising richness in structure, as revealed by
computer simulations. For $K\leq 2$ and for $K\gg 1$, we succeeded in
reproducing this structure as function of $p$ by analytical considerations.

The remainder of this paper is divided into the following parts: In Sec. II we
define the RBN model and the type of noise used for our study. In Sec. III, we
define the first order parameter, the Hamming distance, and evaluate it
numerically and analytically. In Sec. IV, we define the second order
parameter, the frozenness, and evaluate it using computer simulations and
analytical considerations. Finally, we summarize and discuss our findings in
Sec. V.

##  II Model

A Boolean network is defined as a directed network of $N$ nodes representing
Boolean variables $\bm{\sigma}\in\\{1,0\\}^{N}$, which are subject to a
dynamical update rule,

$\bm{\sigma}(t+1)=\bm{f}\left(\bm{\sigma}(t)\right),$ (1)

where $f_{i}$ is a function assigned to node $i$ that depends exclusively on
the states of its inputs.

We introduce noise into the system through a probability $p$ that a node does
not obey its deterministic update rule,

$\bm{\sigma}(t+1)=\bm{f}\left(\bm{\sigma}(t)\right)\veebar\bm{n},$ (2)

where $\bm{n}$ is a random vector, with elements $n_{i}$ being $1$ with
probability $p$ and $0$ otherwise. The symbol $\veebar$ represents the
“exclusive or” Boolean operation. Hence, for $p=0$ the deterministic behaviour
is recovered, and for $p=1/2$ the dynamics is completely stochastic.

RBNs are a special case of Boolean networks, where all possible Boolean
functions are assigned randomly to each node with the same probability, and
where the nodes are randomly connected. The number of inputs of each node is
fixed at a value $K$. The random wiring leads to a Poisson distribution with
mean $K$ for the number of outputs. When updated deterministically, RBNs are
in the frozen phase for $K=1$. After a transient time, they reach an attractor
where all nodes (or all nodes apart from a small number) are permanently
frozen in one of the two Boolean states. Networks with larger $K$ have also a
_frozen core_ of nodes for $p=0$, and the nodes belonging to it become frozen
after a transient time. For $K=2$, all but of the order of $N^{2/3}$ nodes
belong to the frozen core. With increasing $K$, the frozen core contains an
ever smaller proportion of nodes. For $K=2$, the nonfrozen part of the network
consists of several independent components. Each of these components contains
a set of _relevant nodes_ , which are connected such that there is at least
one feedback loop among them, and “trees” of nonfrozen nodes which are rooted
in the relevant nodes and which are slaved to the dynamics of the relevant
nodes.

##  III Average Hamming distance

###  III.1 Definition

We use the average in time of the Hamming distance between the states of two
copies of a network in order to quantify the effect of noise on the dynamics.
Consider a given network in the initial state $\bm{\sigma}(t=0)$, and an exact
replica, which is initially in the same state,
$\bm{\sigma}^{\prime}(t=0)=\bm{\sigma}(t=0)$. The dynamics of both networks
are evolved in parallel, but noise is applied only to
$\bm{\sigma}^{\prime}(t)$, as in Eq. (2). The mean Hamming distance $h(t)$
between the two networks is defined as

$h(t)=\frac{1}{N}\langle\left|\sigma_{i}(t)-\sigma_{i}^{\prime}(t)\right|\rangle\,,$
(3)

where $\langle\dots\rangle$ denotes the average over the noise.

The long-time average $h$ of the Hamming distance is defined as

$h=\lim_{T\to\infty}\frac{1}{T}\frac{1}{N}\sum_{t,i}\left|\sigma_{i}(t)-\sigma_{i}^{\prime}(t)\right|.$
(4)

If the trajectories become completely uncorrelated after some time, we have
$h=1/2$. If the trajectories remain closer in state space, we have $h<1/2$.
The case $h>1/2$ does not occur in our model and is therefore not considered
in this paper.

###  III.2 Annealed approximation

We will first evaluate analytically the Hamming distance by using the so-
called annealed approximation Derrida and Pomeau (1986). This is a mean field
theory, which neglects correlations between nodes and the finite size of the
network. The annealed approximation corresponds to the behaviour of a
(infinitely large) network where all the edges are randomly rewired at each
time step. Within the annealed approximation, the dynamics of a RBN without
noise (i.e. for $p=0$) is fully specified by the parameter $\lambda$, which is
$K$ times the probability that a node changes its state when one (or more) of
its inputs is flipped. For RBNs, we have $\lambda=K/2$, since for any input
combination, there is an equal probability that the output of a function will
be either $0$ or $1$. When considering two replicas of a network, $\lambda$ is
identical to the mean number of nodes that assume a different state in the two
networks at time $t=1$ when at time $t=0$ the state of only one node was
different.

At any time, the Hamming distance between a network with noise and its twin
noiseless counterpart, as described by Eq. (3), is simply the fraction of
nodes which were changed by noise or by the effect of previously changed
nodes. The time evolution of $h(t)$ can then be described as the evolution of
the population of flipped nodes. (A node in the replica with noise is called
“flipped” it its state deviates from the state it has in the replica without
noise.) Let $q(h(t))=1-(1-h(t))^{K}$ denote the probability that a node has at
least one flipped input. Then the probability $h(t+1)$ that a node is flipped
at time $t+1$ can be written as

$\begin{split}h(t+1)&=\frac{\lambda}{K}q(1-p)+\left(1-\frac{\lambda}{K}\right)pq+(1-q)p\\\
&=\frac{\lambda(1-2p)}{K}\left[1-\left(1-h(t)\right)^{K}\right]+p,\\\
\end{split}$ (5)

where the first term in the first line corresponds to the proportion of nodes
that are flipped by previously flipped nodes (and are not flipped back by
noise), and the second and third term are the proportion of nodes that are
flipped by noise (with or without inputs being flipped). The fixed point of
Eq. (5) determines the order parameter $h$, for given $K$ and $p$. We
evaluated this fixed point numerically. Fig. 1 shows $h$ as function of the
noise strength $p$ for several values of $K$. The solid lines are the fixed
point solutions of Eq. (5), the symbols represent the result of computer
simulations of quenched RBNs. The agreement between the annealed approximation
and the real networks is very good.

The most striking feature of Fig. 1 is the existence of a first-order
transition at $p=0$ for $K>2$. This is due to the phase transition to
“chaotic” behaviour for $K>2$. In chaotic networks, even the smallest local
perturbations have a global effect.

###  III.3 The Hamming distance on subsets of nodes

We next evaluate separately the Hamming distance for the frozen core and for
the nonfrozen part of the network. Fig. 2 shows the long-time Hamming
distance, evaluated only for the nodes that belong to the frozen core. These
curves can be fitted using the annealed approximation Eq. (5) under the
condition that the factor $\lambda/K$ on the right-hand side of Eq. (5)
(representing the probability that a node is flipped when at least one input
is flipped) is replaced with $\lambda_{\text{eff}}/K$, with
$\lambda_{\text{eff}}$ being used as a fit parameter. For $K=1$ and $2$, the
frozen core is virtually indistinguishable from the rest of the network, and
$\lambda_{\text{eff}}=\lambda$, but for $K>2$, $\lambda_{\text{eff}}$
decreases with increasing $K$. The reason is that the frozen core becomes
composed mainly of nodes with constant functions which have
$\lambda_{\text{eff}}=0$.

Figure 1: (Color online) Average Hamming distance as function of noise
strength for RBNs of size $N=10^{4}$ for different values of $K$. Each point
was obtained by averaging the results over 3 different network realizations.
The solid lines are the respective steady-state solutions of Eq. (5).

Before evaluating $h$ for the nonfrozen nodes, let us consider the simplest
possible connected set of nonfrozen nodes, which is a simple loop. For $K=1$
and $K=2$, such simple loops of nonfrozen nodes play an important role at
determining the attractors with periods larger than 1 Kaufman and Drossel
(2005), however, a considerable fraction of $K=2$ networks also have more
complex relevant components. The effect of noise on such loops is very
different from its effect on the frozen core, since if one of its nodes is
flipped, this flip propagates indefinitely around the loop. One can calculate
the accumulation of flips on such loops by considering the average Hamming
distance at a given time between a loop without and with noise,

$\begin{split}h(t)&=\sum_{i=0}^{\lfloor\frac{t-1}{2}\rfloor}{t\choose
2i+1}p^{2i+1}(1-p)^{t-(2i+1)}\\\
&\approx\sum_{i=0}^{\lfloor\frac{t-1}{2}\rfloor}\frac{(tp)^{2i+1}e^{-tp}}{(2i+1)!}\\\
&\approx\frac{1}{2}\left(1-e^{-2tp}\right)\,.\end{split}$ (6)

The first equation evaluates the probability that a node has been flipped an
odd number of times, and the subsequent transformations are valid for $t\gg
1$. The Hamming distance approaches the value $1/2$ with an exponential decay,
and with an characteristic time $\tau=1/(2p)$.

Figure 2: (Color online) Average Hamming distance as function of noise for RBN
of size $N=10^{4}$ for different values of $K$. Only nodes belonging to the
frozen core of the network (without noise) were considered. Each point on each
curve was obtained by averaging the results for 3 different network
realizations. The solid lines are the respective solutions of Eq. (5), with
$\lambda_{\text{eff}}$ being used as a fit parameter.

We evaluated how fast a trajectory leaves an attractor in the presence of
noise by first letting the system approach an attractor and by then turning on
the noise and measuring the Hamming distance $h(l)$ to the initial state after
one attractor period $l$. This is identical to the distance from the state of
the noiseless replica, which returns to the initial state at time $l$. Fig. 3
shows the values of $h(l)$ for RBNs with different values of $K$. For $K=1$,
the data match Eq. (6) very well, since the nonfrozen part of the network in
this case can only be composed of simple loops. For $K\geq 2$, the data points
are considerably above this exponential curve because a node can become
flipped via many different paths. The data are better fitted using Eq. (5), in
particular for long periods (i.e. large times). Just as for the case of the
frozen core, $\lambda_{\text{eff}}$ was used as a fit parameter.

For smaller values of the attractor period, the data are considerably below
the fitted line. The reason is that these attractor periods are much smaller
than typical attractor periods, and networks with such short attractors are
not characteristic of the ensemble, but have a state-space structure with a
smaller set of recurrent states. Consequently trajectories diverge less fast
than in typical networks.

Figure 3: (Color online) Average Hamming distance $h(l)$ of the relevant
components of RBNs with different values of $K$, after a full period $l$, with
$p=0.01$. The curves were obtained by sampling at least $2\times 10^{4}$
attractors of several distinct RBNs of sizes $N=10^{2}$, $50$ and $25$ for
$K\leq 2$, $3$ and $4$, respectively. The solid lines are given by Eq. 6 for
$K>1$ and Eq. 5 for $K\geq 2$.

##  IV Frozenness

###  IV.1 Definition

The “frozenness” of a network measures the extent to which the nodes spend
more time in one of the two Boolean states. It is zero, when the nodes spend
the same time in both states, and it is 1 when the network is frozen. The
frozenness of node $i$ is defined by the expression

$\Omega_{i}=\left(q_{0}^{(i)}-q_{1}^{(i)}\right)^{2},$ (7)

where $q_{\sigma}^{(i)}$ is the proportion of time node $i$ is in state
$\sigma$,

$q^{(i)}_{\sigma}=\lim_{T\to\infty}\frac{1}{T}\sum_{t=0}^{T}\delta_{\sigma_{i}(t),\sigma}.$
(8)

By eliminating one of the two probabilities from Eq. (7), we obtain

$\Omega_{i}=\left(2q_{\sigma}^{(i)}-1\right)^{2},$ (9)

where $\sigma$ is either $0$ or $1$.

The frozenness $\langle\Omega\rangle$ of the network is obtained by averaging
over the nodes. Figure 4 shows the frozenness $\langle\Omega\rangle$ as
function of the noise strength $p$ obtained by computer simulations of
networks of size $10^{4}$, for different $K$.

Figure 4: (Color online) Frozenness as function of noise strength for RBNs of
size $N=10^{4}$ for different values of $K$. Each curve was obtained by
averaging the results for 100 different network realizations. The solid lines
correspond to the averages of $\rho_{\Omega}(\Omega|K,p)$, obtained in Secs.
IV.2, IV.3 and IV.4, for $K=1$, $2$ and $\geq 3$, respectively.

At $p=0$, the frozenness corresponds obviously to the size of the frozen core,
and it decreases towards 0 as the noise strength approaches the value 0.5.

In order to derive these curves analytically, the annealed approximation is of
no use, since a network that is rewired during the course of time has very
small frozenness, which is due uniquely to the constant functions. Therefore,
a simple analytical calculation, which does not require the consideration of
correlations between nodes, can only be performed for nodes with constant
functions: for such nodes the frozenness is given by

$\Omega_{i}=(1-2p)^{2}\,.$ (10)

In the following, we will present more advanced analytical evaluations and
further computer simulations for RBNs with different values of $K$.

###  IV.2 $K=1$

For $K=1$ there are only $2^{2^{1}}=4$ possible Boolean functions, two of
which are constant ($1$ or $0$), and the remaining ones are the copy
($f(\sigma)=\sigma$) and invert ($f(\sigma)=\lnot\sigma$) functions. As far as
the analysis of frozenness is concerned, there are only two distinct
functions, constant and non-constant, since the output value is not relevant,
but only how often it changes. Since each of the two types of functions occurs
equally often in a $K=1$ RBN, we have $\lambda=1/2$, but $K=1$ networks with
other values of $\lambda$ can also be constructed.

In a network with $K=1$, each node has only one input, and this input node
also has one input, etc. In order to evaluate the probability that a node is
flipped, one only needs to consider the chain of those nodes that can have an
influence on the considered node. Nodes with constant functions present a
barrier to the propagation of a perturbation, since they do not respond to a
change in their inputs, and therefore the chain ends (or, more precisely,
begins) at a node with a constant function.

Without loss of generality, we define $q^{(i)}$ as being the proportion of
time node $i$ assumes its most frequent value,

$q^{(i)}\equiv\max(q^{(i)}_{\sigma},1-q^{(i)}_{\sigma})\;\in\left[1/2,1\right].$
(11)

Since the value of $q^{(i)}$ is fully determined by the distance of node $i$
to a node with a constant function, we choose the label $i$ in the remainder
of this subsection to signify this distance. If the node itself has a constant
function, we have $i=0$, if the node has a non-constant function, but its
input has a constant function, we have $i=1$, etc.

The value of $q^{(0)}$, i.e. for nodes with constant functions, is simply

$q^{(0)}(p)=1-p.$ (12)

For larger values of $i$, we have the recursion relation

$\begin{split}q^{(i)}&=(1-p)q^{(i-1)}+(1-q^{(i-1)})p\\\
&=\frac{1}{2}(1-2p)^{i+1}+\frac{1}{2},\end{split}$ (13)

where the solution of the recursion relation was obtained using Eq. (12). The
probability of finding a given $q^{(i)}$ in the network is

$p_{q}(q^{(i)}|p)=(1-\lambda)\lambda^{i}=\left(\frac{1}{2}\right)^{i+1}\,.$
(14)

The frozenness of the network is thus given by

$\left<\Omega\right>(p)=\sum_{i}p_{q}(q^{(i)}|p)(1-2q^{(i)})^{2}\,,$ (15)

which is plotted in Figure 4 and fits the curve for $K=1$ well.

Although networks with $K=1$ have only a discrete set of possible $q$-values,
the distribution of $q$ values appears as a continuum when determined by
computer simulations. There are two reasons for this. First, the data points
for $q$ close to 1/2 (i.e. for $\Omega$ close to 0) are so close to each other
that they cannot be resolved, since a computer simulation uses a non-vanishing
bin size. Solving Eq. (13) for $i$, inserting the result in (14), and using
the relation $\rho_{q}(q|K=1,p)dq=p_{q}(q|K=1,p)$ with $dq=q^{(i)}-q^{(i+1)}$,
we obtain

$\rho_{q}(q|K=1,p)\propto(2q-1)^{\frac{\ln\lambda}{\ln(1-2p)}-1}\,,$ (16)

with $\lambda=1/2$ for RBNs. Thus, in the limit $q\to 1/2$ the distribution of
$2q-1$ follows as a power-law with an exponent given by the above expression.
The probability density of $\Omega$ also decays as a power law, in the limit
$\Omega\to 0$, but with a different exponent, since

$\begin{split}\rho_{\Omega}(\Omega|K=1,p)=&\;\rho_{q}\left(\left.\frac{\sqrt{\Omega}+1}{2}\right|K=1,p\right)\frac{dq}{d\Omega}\\\
\sim&\;\Omega^{\frac{\ln\lambda}{\ln(1-2p)}-3/2}\,.\end{split}$ (17)

Second, a computer simulation averages only over a finite amount of time, $T$,
and therefore the measured values $q^{\prime}$ are Gaussian distributed around
the exact value $q$,

$\begin{split}\tilde{\rho}_{q^{\prime}}(q^{\prime}|K=1,T,q)\approx\frac{1}{\sqrt{2\pi
q/T}}e^{-\frac{(q^{\prime}-q)^{2}}{2q/T}},\end{split}$ (18)

This dependence on $T$ can be included in Eq. (14) to obtain the probability
density function for $q$

$\begin{split}\rho_{q}(q|K=1,p,T)=&\,(1-\lambda)\sum_{i=0}^{\infty}\lambda^{i}\tilde{\rho}_{q^{\prime}}(q|T,q^{(i)}).\end{split}$
(19)

For the distribution of $\Omega$ values, we obtain

$\rho_{\Omega}(\Omega|K,p,T)=\rho_{q}\left(\left.\frac{\sqrt{\Omega}+1}{2}\right|K,p,T\right)\frac{1}{2\sqrt{\Omega}}.$
(20)

Fig. 5 shows the distribution of the frozenness for a quenched network with
$K=1$, for several values of $p$. It can be seen that there is very good
agreement with Eq. (20). The presence of fluctuations significantly deviates
some of the distributions from the expected power-law decay. For $p=0.01$, the
small values of frozenness, which are not in agreement with the theoretical
result, are due to the existence of loops, which are omitted in the analysis
above. The probability that a node is part of a nonfrozen loop tends to zero
as the network becomes larger, and therefore these points vanish in the limit
of infinite system size.

Figure 5: (Color online) Distribution of the frozenness for different values
of the noise strength for RBNs of size $N=10^{5}$ for $K=1$ and $T=10^{4}$.
Each curve was obtained by averaging over 100 different network realizations.
The solid lines are given by Eq. (20).

###  IV.3 $K=2$

As $K$ becomes larger, the number of possible functions grows very fast as
$2^{2^{K}}$, and a detailed analysis of the frozenness, as was done for $K=1$,
becomes more complicated. The values of $q$ are still discontinuously
distributed, but their number increases fast with $K$, due to the numerous
combinations of Boolean functions that can determine the $q$-values of the $K$
inputs of a node and thus, in combination with the node’s Boolean function,
the $q$-value of its output. Here we will lay out the basic considerations
needed to obtain the distribution of $q$ for all $K>1$, and we will obtain by
numerical iteration the distribution for $K=2$. Without loss of generality, we
redefine $q$ as $q\equiv q_{\sigma=1}$, i.e., the fraction of time a given
node has the value 1 (as opposed to Eq. (11), which simplified the case
$K=1$).

In general, the probability density function $\rho_{q}(q|K=1,p)$ needs to
account for all possible recursive combinations of output functions and their
inputs. We can thus write the following self-consistent expression,

$\rho_{q}(q|K,p)=\int_{0}^{1}\dotsi\int_{0}^{1}\sum_{f}p_{f}\delta(q^{(f)}-q)\prod_{i=1}^{K}\rho_{q}(q^{(i)}|K,p)dq^{(i)},$
(21)

where the sum is taken over all Boolean functions; $p_{f}$ is the probability
of the $f$th Boolean function ($p_{f}=2^{-2^{k}}$ ), and

$q^{(f)}=(1-2p)q^{(f)}(\\{q^{(i)}\\})+p,$ (22)

where $q^{(f)}(\\{q^{(i)}\\})$ is the value of $q$ for a specific function
$f$, given the values $\\{q^{(i)}\\}$ of its inputs, for $i=1,\dots,K$.

Since Eq. (21) involves an expression $q^{(f)}(\\{q^{(i)}\\})$ for all Boolean
functions, a general closed solution becomes unfeasible. However, for $K=2$
Eq. (21) can at least be solved numerically, since there are only $16$
possible functions, given in Table 1. Eq. (21) is then solved by iteration,
until convergence to a self-consistent $q$ distribution is obtained. We
started with the initial distribution

$\rho^{0}_{q}(q|K,p)=\frac{1}{2}\delta(q-p)+\frac{1}{2}\delta(q-(1-p))\,.$
(23)

In the end, we determined the final distribution $\rho_{\Omega}(\Omega|K=2,p)$
by using Eq. (20).

$f$ | $f_{i}(\sigma_{1},\sigma_{2})$ | $q^{(f)}(q_{1},q_{2})$
---|---|---
$0$ | $0$ | $0$
$1$ | $\sigma_{1}\land\sigma_{2}$ | $q_{1}q_{2}$
$2$ | $\sigma_{1}\land\lnot\sigma_{2}$ | $q_{1}(1-q_{2})$
$3$ | $\sigma_{1}$ | $q_{1}$
$4$ | $\lnot\sigma_{1}\land\sigma_{2}$ | $(1-q_{1})q_{2}$
$5$ | $\sigma_{2}$ | $q_{2}$
$6$ | $\sigma_{1}\veebar\sigma_{2}$ | $q_{1}+q_{2}-2q_{1}q_{2}$
$7$ | $\sigma_{1}\lor\sigma_{2}$ | $q_{1}+q_{2}-q_{1}q_{2}$
$8$ | $\lnot(\sigma_{1}\lor\sigma_{2})$ | $1-q^{(7)}(q_{1},q_{2})$
$9$ | $\lnot(\sigma_{1}\veebar\sigma_{2})$ | $1-q^{(6)}(q_{1},q_{2})$
$10$ | $\lnot\sigma_{2}$ | $1-q_{2}$
$11$ | $\sigma_{1}\lor\lnot\sigma_{2}$ | $q^{(7)}(q_{1},1-q_{2})$
$12$ | $\lnot\sigma_{1}$ | $1-q_{1}$
$13$ | $\lnot\sigma_{1}\lor\sigma_{2}$ | $q^{(7)}(1-q_{1},q_{2})$
$14$ | $\lnot(\sigma_{1}\land\sigma_{2})$ | $1-q_{1}q_{2}$
$15$ | $1$ | $1$
Table 1: Expressions of $q^{(f)}(q_{1},q_{2})$ for all Boolean functions for
$K=2$. The Boolean expressions of each function is also given for reference.

Fig. 6 shows the distribution of $\Omega$ for simulated quenched RBNs with
$K=2$ for different values of $p$, compared with the result of the numerical
evaluation of Eq. (21) as described above. There is a very good agreement
between the two types of results. The peaks correspond to prominent values of
the frozenness. The rightmost peak is always due to the constant functions,
but large frozenness values are also obtained for other functions. For
instance, $f_{1}$ assumes the value 0 whenever both inputs are different. If
both inputs have $q=1/2$ (i.e., $\Omega=0$), the value of $q_{f=1}$ is
$(1-2p)1/4+p$ and $\Omega=((1-2p)/2+2p-1)^{2}$. This is the second main peak
of $\rho_{\Omega}(\Omega|2,0.4,10^{4})$ (counted from the right end). For
smaller values of $p$, the peaks are not discernible, and a broad continuum
appears, with a distribution that follows a power-law with an exponent $\simeq
0.6$ as $\Omega\to 0$. As $T$ becomes larger, it is expected that the
continuous regions become more and more discontinuous, as can be seen in Fig.
7, which shows the theoretical prediction for larger times $T$. Moreover, when
the resolution is increased (see inset), it can be seen that peak-like regions
which appear like fluctuations around a single value of $\Omega$, are in fact
composed of sharper peaks, which themselves are composed of other peaks,
building a fractal structure.

Another distinguishing feature seen in Fig. 6 is a sharp transition at $p=0$,
where the only two possible values of $q$ are $0$ and $1$, both of which
amount to $\Omega=1$, leading to the variance $\sigma_{\Omega}^{2}=0$. For
$p>0$ this abruptly changes, and a wide range of values of $q$ are possible,
which discontinuously leads to $\sigma_{\Omega}^{2}>0$. There is no such
discontinuous transition for other $K$ values, but a continuous one (see
following section), which makes the case $K=2$ special.

Figure 6: (Color online) Distribution of the frozenness of nodes for different
noise strength for RBNs of size $N=10^{5}$ and for $K=2$ and $T=10^{4}$. Each
curve was obtained by averaging over 100 different network realizations. The
solid lines show the values of $\rho_{\Omega}(\Omega|K=2,p,10^{4})$ according
to Eq. (20). The line segment corresponds to a power law with an exponent
$0.6$. Figure 7: (Color online) Expected distribution of frozenness of each
node as function of noise for RBNs with $K=2$ and different values of $T$. The
curves represent Eq. (20). The inset shows a zoom into
$\rho_{\Omega}(\Omega|K=2,0.4,T)$ for different values of $T$.

The function $\langle\Omega\rangle$ is obtained by performing the integral
$\int\Omega\cdot\rho_{\Omega}(\Omega|K=2,p,T)d\Omega$. As can be seen in
Figure 4, our calculation of this function agrees well with the results of
computer simulations.

###  IV.4 $K>2$

For larger values of $K$, numerical solutions of Eq. (21) become progressively
more elaborate. We did not pursue the task of writing the expressions of 256
function $q_{f}(\\{q^{(i)}\\})$ for $K=3$ or of $65,536$ functions for $K=4$.
Instead, we perform in the following an approximation that is good for a large
number of inputs per node.

When $K$ is large, the vast majority of Boolean functions have the output $1$
for approximately half the input combinations. This means that almost all
nodes have at their inputs $q$ values close to $1/2$. We therefore make the
assumption that the input values to each node are $1$ and $0$ with probability
$1/2$, independently from each other. This means that for any given function
all input combinations are equally probable. It then follows immediately that
the possible $q$ values are identical to the possible fractions of output
values 1 in the truth table of a Boolean function, and that the probability
for a given $q$ value is

$p_{q}(q|K,p=0)=2^{-M}{M\choose qM}\,.$ (24)

Here, we have defined $M=2^{K}$, and the possible $q$ values are thus
multiples of $1/M$.

In the presence of noise, each output value is inverted with probability $p$,
implying that $q$ is changed to $q^{\prime}=q(1-p)+p(1-q)=(1-2p)q+p$. We
therefore have

$\rho_{q}(q|K,p)=\frac{M}{1-2p}p_{q}\left(\left.\frac{q-p}{1-2p}\right|K,p=0\right)\,.$
(25)

For the frozenness $\Omega$ we obtain the distribution

$\rho_{\Omega}(\Omega|K,p)=\frac{M}{2(1-2p)\sqrt{\Omega}}p_{q}\left(\left.\frac{\sqrt{\Omega}+1-p}{2(1-2p)}\right|K,p\right)\,.$
(26)

Finally, one needs to take into account the effect of fluctuations, exactly as
was done for the previous cases,

$\rho_{q}(q|p,K,T)=\int_{0}^{1}\rho_{q}(q^{\prime}|p,K)\tilde{\rho}_{q}(q|T,q^{\prime})dq^{\prime},$
(27)

where $\tilde{\rho}_{q}(q|T,q^{\prime})$ is given by Eq. (18).

Fig. 8 shows the distributions of frozenness for $K=4$ and $5$. In contrast to
the cases $K=1$ and $K=2$, the peaks are less pronounced, and are hardly
visible. For $K=3$ (not shown) there are some peaks which are still visible,
specially for high values of $p$. Therefore, the high-$K$ approximation Eq
(27) is very good already for $K=4$. Both distributions show the same power-
law decay $\rho_{\Omega}(\Omega|K,p)\sim\Omega^{-1/2}$. This is simply due to
the fact that for $\Omega\to 0$ ($q\to 1/2$) the shape of $\rho_{q}(q|K,p)$ is
essentially flat, and thus $\rho_{\Omega}(\Omega|K,p)\sim
dq/d\Omega=\Omega^{-1/2}/2$.

Figure 8: (Color online) Distribution of the frozenness of nodes as function
of noise strength for RBN of size $N=10^{4}$ for different values of $K$. Each
curve was obtained by averaging the results for 100 different network
realizations. The solid lines correspond to Eq. 27.

The quality of our approximation can also be assessed by comparing the
analytical prediction for $\left<\Omega\right>$ with the computer simulations
(Fig. 4). We expect that the approximation becomes even better for larger $K$.

##  V Conclusion

We have investigated the effect of thermal noise on RBNs by evaluating two
order parameters, the long time average of the Hamming distance between two
networks and the average frozenness of the network. While for $K=1$ and $K=2$
the average Hamming distance increases continuously from 0 to 1/2 as $p$
increases from 0 to 1/2, it has a jump at $p=0$ for $K>2$. These findings are
well reproduced by the annealed approximation, and they are a consequence of
the transition from a frozen to a “chaotic” phase in the deterministic system.
In the chaotic phase (occuring for $K>2$), initially nearby trajectories
become eventually uncorrelated. The smooth increase of the Hamming distance
towards the value 1/2 is compatible with what was found in Miranda and Parga
(1989); Qu et al. (2002) for small networks.

The analysis of the average frozenness of the network required more
sophisticated calculations than the annealed approximation, and revealed
intricate details of the network dynamics. For all values of $K$ the
probability distribution of frozenness is a set of delta peaks. For $K=1$,
these peaks can be obtained by considering the distance of nodes to nodes with
constant functions. For $K>1$, the analysis becomes a lot more elaborate, due
to the large number of Boolean functions and the resulting vast number of
possible combinations of frozenness values for the inputs of each function. We
explained the general method, and performed the actual numerical evaluation
for the case $K=2$. The delta peaks show a fractal structure, which emerges
from the iterated recursion relation for the possible frozenness values. The
variance of the frozenness distribution changes continuously with $p$ for all
$K\neq 2$, but for $K=2$ it has a jump at $p=0$, where the variance changes
discontinuously from $0$ to a value larger than zero. For larger values of
$K$, the delta peaks are so close to each other that the frozenness
distribution appears continuous, and in this limit we succeeded in performing
an approximate analytical calculation.

We do not find a phase transition at finite noise strength, in contrast to
Huepe and Aldana-González (2002), where Boolean networks with threshold
functions following a majority rule were used. Such a system undergoes a
second order phase transition from an ordered “ferromagnetic” phase, where all
nodes assume the same value for the majority of time, to a disordered phase,
where the nodes assume both states equally often. The presence of an ordered
phase is a direct consequence of the majority rule, and this transition is
similar to that in a network of Ising spins Indekeu (2004); Aleksiejuk et al.
(2002). The order parameter in Huepe and Aldana-González (2002) was defined as
the average “alignment” $s=\left|\left<1-2\sigma\right>\right|$, which is $1$
if all nodes are in the same state. The order parameter $s$ is only meaningful
in systems where the system is ordered in the absence of noise, and where the
symmetry between the states with values $0$ and $1$ is broken, as in
ferromagnetic spin systems. Otherwise, $\left<\Omega\right>$ is a better order
parameter, because it captures disordered frozen phases, such as for $K=1$ in
RBNs. Of course, a phase transition in the value of $s$ is always accompanied
by a phase transition in the value of $\left<\Omega\right>$. The opposite is
not always true.

It is to be expected that real networks show some kind of robustness to noise,
since they must be able to carry out their function in a noisy environment. As
the results of this work show, only for RBNs with $K=1$ do the order
parameters change slowly as noise is switched on. RBNs with $K>1$ fail to
exhibit robustness to noise, which is hardly surprising given the random
wiring of the system and the random choice of functions. It will therefore be
interesting to extend the present study to networks with a more restricted set
of functions with more biological relevance, such as threshold Szejka et al.
(2008) or canalizing functions Moreira and Amaral (2005); Kauffman et al.
(2004). At least for some sets of functions, one should expect a phase
transition at a finite noise strength, similar to the transition seen in Huepe
and Aldana-González (2002). The survival of the “ordered” phase up to a
certain noise strength can be viewed as a certain type of robustness.

It remains to be seen how other network topologies Iguchi et al. (2007) and
the incorporation of redundancy Gershenson et al. (2006) change a network’s
response to noise. In Iguchi et al. (2007), it was shown that a scale-free
input distribution changes the average number and length of attractors. In
Gershenson et al. (2006), redundancy was introduced as functional duplications
of nodes in the network, which resulted in greater robustness against random
mutations of the update functions. In both papers, only deterministic dynamics
were considered. The effects of these (or other more general) topological and
functional characteristics may strongly alter the response of a network to
thermal noise. Finding the general conditions required for reliable dynamics
in a stochastic environment will be an important step towards a deeper
understanding of the dynamical features of real networks.

We acknowledge the support of this work by the Humboldt Foundation.

## References

  * Kauffman (1969) S. A. Kauffman, J. Theor. Biol. 22, 437 (1969).
  * Drossel (2008) B. Drossel, _Reviews of Nonlinear Dynamics and Complexity_ (Wiley, 2008), vol. 1, ISBN 3527407294.
  * Rosen-Zvi et al. (2001) M. Rosen-Zvi, A. Engel, and I. Kanter, Phys. Rev. Lett. 87, 078101 (2001).
  * Moreira et al. (2004) A. A. Moreira, A. Mathur, D. Diermeier, and L. A. N. Amaral, Proc. Nat. Ac. Sci. 101, 12085 (2004).
  * Lagomarsino et al. (2005) M. C. Lagomarsino, P. Jona, and B. Bassetti, Phys. Rev. Lett. 95, 158701 (2005).
  * Bornholdt (2005) S. Bornholdt, Science 310, 449 (2005).
  * Greil and Drossel (2005) F. Greil and B. Drossel, Phys. Rev. Lett. 95, 048701 (2005).
  * Klemm and Bornholdt (2005) K. Klemm and S. Bornholdt, Phys. Rev. E 72, 055101 (2005).
  * Samuelsson and Troein (2003) B. Samuelsson and C. Troein, Phys. Rev. Lett. 90, 098701 (2003).
  * Drossel (2005) B. Drossel, Phys. Rev. E 72, 016110 (2005).
  * Drossel et al. (2005) B. Drossel, T. Mihaljev, and F. Greil, Phys. Rev. Lett. 94, 088701 (2005).
  * Shmulevich et al. (2002) I. Shmulevich, E. R. Dougherty, S. Kim, and W. Zhang, Bioinformatics 18, 261 (2002).
  * Indekeu (2004) J. O. Indekeu, Physica A 333, 461 (2004).
  * Aleksiejuk et al. (2002) A. Aleksiejuk, J. A. Holyst, and D. Stauffer, Physica A 310, 260 (2002).
  * McAdams and Arkin (1997) H. H. McAdams and A. Arkin, Proc. Nat. Ac. Sci. 94, 814 (1997).
  * Fretter and Drossel (2008) C. Fretter and B. Drossel, Eur. Phys. J. B 62, 365 (2008).
  * Miranda and Parga (1989) E. N. Miranda and N. Parga, Europhys. Lett. 10, 293 (1989).
  * Golinelli and Derrida (1989) O. Golinelli and B. Derrida, J. Phys 50, 1587 (1989).
  * Qu et al. (2002) X. Qu, M. Aldana, and L. P. Kadanoff, J. Stat. Phys. 109, 967 (2002).
  * Huepe and Aldana-González (2002) Huepe and Aldana-González, J. Stat. Phys. 108, 527 (2002).
  * Derrida and Pomeau (1986) B. Derrida and Y. Pomeau, Europhys. Lett. 1, 45 (1986), ISSN 0295-5075.
  * Kaufman and Drossel (2005) Kaufman and Drossel, Eur. Phys. J. B 43, 115 (2005).
  * Szejka et al. (2008) A. Szejka, T. Mihaljev, and B. Drossel, New J. Phys. 10, 063009 (2008).
  * Moreira and Amaral (2005) A. A. Moreira and L. A. N. Amaral, Phys. Rev. Lett. 94, 218702 (2005).
  * Kauffman et al. (2004) S. Kauffman, C. Peterson, B. Samuelsson, and C. Troein, Proc. Nat. Ac. Sci. 101, 17102 (2004).
  * Iguchi et al. (2007) K. Iguchi, S. ichi Kinoshita, and H. S. Yamada, J. Theor. Biol. 247, 138 (2007).
  * Gershenson et al. (2006) C. Gershenson, S. A. Kauffman, and I. Shmulevich, _Artificial Life X: Proceedings of the Tenth International Conference on the Simulation and Synthesis of Living Systems_ (The MIT Press, 2006), ISBN 0262681625\.


# A Quantum Finite Automata Approach to Modeling the Chemical Reactions

Amandeep Singh Bhatia${}^{1,^{*}}$, Shenggen Zheng2
1Chitkara University Institute of Engineering & Technology, Chitkara
University, Punjab, India
2Center for Quantum Computing, Peng Cheng Laboratory, Shenzhen, China
E-mail<EMAIL_ADDRESS>
###### Abstract

In recent years, the modeling interest has increased significantly from the
molecular level to the atomic and quantum scale. The field of computational
chemistry plays a significant role in designing computational models for the
operation and simulation of systems ranging from atoms and molecules to
industrial-scale processes. It is influenced by a tremendous increase in
computing power and the efficiency of algorithms. The representation of
chemical reactions using classical automata theory in thermodynamic terms had
a great influence on computer science. The study of chemical information
processing with quantum computational models is a natural goal. In this paper,
we have modeled chemical reactions using two-way quantum finite automata,
which are halted in linear time. Additionally, classical pushdown automata can
be designed for such chemical reactions with multiple stacks. It has been
proven that computational versatility can be increased by combining chemical
accept/reject signatures and quantum automata models.

Keywords: chemical reaction, quantum finite automata, pushdown automata, two-
way quantum finite automata, Belousov-Zhabotinsky reaction

##  I Introduction & Motivation

Recently, the connection between complex reactions and their thermodynamics
has got overwhelming response among research communities. Initially, in the
1970s, Conrad (1, ) processed the information of molecular systems and stated
that complex biochemical systems cannot be analyzed on classical computers.
Till now, artificial approaches use complex biomolecules or logic gates based
reaction-diffusion systems to solve the problems 23 ; 24 ; 25 . Classical
systems are not robust and capable to describe quantum systems. Some tasks
that are impossible in classical systems can be realized in quantum systems.
Quantum computing is concerned with computer technology based on the
principles of quantum mechanics, which describes the behavior and nature of
matter and energy on the quantum level (2, ). Quantum computation demonstrates
the computation power and other properties of the computers based on
principles of quantum mechanics.

Models of finite automata are abstract computing devices, which play a crucial
role to solve computational problems in theoretical computer science.
Classical automata theory is closely associated with formal language theory,
where automata are ranked from simplest to most powerful based on their
language recognition power 3 . Classical automata theory has been of
significant importance due to its practical real-time applications in the
development of several fields. Therefore, it is a natural goal to study
quantum variants of classical automata models, which play an important role in
quantum information processing.

The quantum automata theory has been developed using the principles of quantum
mechanics and classical automata. Quantum computational models make it
possible to examine the resources needed for computations. Soon after the
brainstorm of Shor’s factorization quantum algorithm 4 , the first models of
quantum finite automata (QFAs) have been introduced. Initially, Kondacs and
Watrous 6 , and Moore and Crutchfield 5  proposed the concept of quantum
automata separately. There is a diversity of quantum automata models have been
studied since then and demonstrated in various directions such as quantum
finite automata 5 , Latvian QFA, 1.5-way QFA, two-way QFA (2QFA) 6 , quantum
finite state machines of matrix product states 54 , quantum sequential
machines 58 , quantum pushdown automata, quantum Turing machine, quantum queue
automata 56 , quantum multicounter machines (QMCM), quantum multihead finite
automata (QMFA), quantum finite automata with classical states (2QCFA) 36 ; 37
, quantum omega automata 55 , state succinctness of two-way probabilistic
finite automata (2PFA), QFA, 2QFA and 2QCFA 51 ; 48 ; 49 , interactive proof
systems with quantum finite automata 52 ; 53  and semi-quantum two-way finite
automata 46 ; 47  and many more since last two decades (50, ; 57, ; 8, ; 9, ;
10, ). These models are effective in determining the boundaries of various
computational features and expressive power. Quantum computers are more
powerful than Turing machines and even probabilistic Turing machines. Thus,
mathematical models of quantum computation can be view as generalizations of
its physical models.

Computational biochemistry has been a rapidly evolving research area at the
interface between biology, chemistry, computer science, and mathematics. It
helps us to apply computational models to understand biochemical and chemical
processes and their properties. A combination of chemistry and classical
automata theory provides a constructive means of refining the number of
objects that allow to understand the energetic cost of computation 11 . The
research has been consistently grown in the field of chemical computing. There
exist two ways to model complex chemical reactions: abstract devices and
formal models based on multiset rewriting 13 . The complex chemical reaction
networks carry out chemical processes that mimic the workings of classical
automata models. Recently, Duenas-Diez and Perez-Mercader 11 ; 12  have
designed chemical finite automata for regular languages and chemical automata
with multiple stacks for context-free and context-sensitive languages.
Further, the thermodynamic interpretation of the acceptance/rejection of
chemical automata is given. It is useful to understand the energetic cost of
chemical computation. They have used the one-pot reactor (mixed container),
where chemical reactions and molecular recognition takes place after several
steps, without utilizing any auxiliary geometrical aid.

In classical automata theory, it is known that two-way deterministic finite
automata (2DFA) can be designed for all regular languages. It has also been
investigated that two-way probabilistic finite automata (2PFA) can be designed
for a non-regular language $L=\\{a^{n}b^{n}\mid n\geq 1\\}$ in an exponential
time 14 ; 15 . The research has consistently evolved in the field of quantum
computation and information processing. In quantum automata theory, it has
been proved that 2QFA can be designed for L with one-sided bounded error and
halted in linear time. Moreover, it has been demonstrated that 2QFA can be
also designed for non-context free language $L=\\{a^{n}b^{n}c^{n}\mid n\geq
1\\}$ 6 . Hence, 2QFA is strictly more powerful than its classical
counterparts based on language recognition capability.

The field of chemistry and chemical computing plays a significant role in the
evolution of computational models to mimic the behavior of systems at the
atomic level. It is greatly influenced by the computing power of quantum
computers. Motivated from the above-mentioned facts, we have modeled chemical
reactions in the form of formal languages and represented them using two-way
quantum finite automata. The main objective is to examine how chemical
reactions perform chemical sequence identification equivalent to quantum
automata models without involving bio-chemistry or any auxiliary device. The
crucial advantage of this approach is that chemical reactions in the form of
accept/reject signatures can be processed in linear time with one-sided
bounded error. The organization of rest of this paper is as follows:
Subsection is devoted to prior work. In Sect. 2, some preliminaries are given.
The definition of two-way quantum finite automata is given in Sect. 3. In
Sect. 4, the chemical reactions are transcribed in formal languages and
modeled using two-way quantum finite automata approach. Summary of work is
given in Section 5. Finally, Sect. 6 is the conclusion.

###  I.1 Prior work

The field of chemical computation has rich and interesting history. Various
researchers have represented the chemical computation using the concept of
logic gates based reaction-diffusion systems and artificial intelligence
approaches. In early 1970s, Conrad 1  differentiated the information
processing on molecules from digital computing. Nearly a decade after, Okamoto
et al. 16  proposed the concept of a theoretical chemical diode for cyclic
enzyme systems. It has been proved that it can be used to analyze the dynamic
behavior of metabolic switching events in bio-computer. In 1991, Hjelmfelt et
al. 17  designed neural networks and finite state machines using chemical
diodes. It has been investigated that the execution of a universal Turing
machine is possible using connecting chemical diodes.

In 1995, Tóth and Showalter 18  implemented AND and OR logic gates using
reaction-diffusion systems, where the signals are programmed by chemical
waves. It was the first empirical realization of chemical logical gates. In
1997, Magnasco 19  showed that logic gates can be constructed and executed in
the chemical kinetics of homogeneous solutions. It has been proved that such
constructions have computational power equivalent to Turing machine. Adamatzky
and Lacy Costello 20  experimentally realized the Chemical XOR gate by
following the same approach of Toth and Showlter in 2002. Further, Górecki et
al. 21  constructed the chemical counters for an information processing in the
excitable reaction-diffusion systems.

It is one of the most promising new areas of research. Some of the
difficulties can be caused by connecting several gates together for advanced
computation. Thus, recently, researchers started focusing on native chemical
computation, i.e. without reaction-diffusion systems. In 1994, Adleman 22
proposed the concept of DNA computing and solved the Hamiltonian path problem
by changing DNA strands. In 2009, Benenson 23  reviewed biological measurement
tools for new generation biocomputers. Prohaska et al. 24  studied protein
domain using chromatin computation and introduced chromatin as a powerful
machine for chemical computation and information processing. In 2012, Bryant
25  proved chromatin computer as computationally universal and used to solve
an example of combinatorial problem.

The structures of DNA and RNA are represented using the concept of classical
automata theory 26  27 . Krasinski et al. 28  represented the restricted
enzyme in DNA with pushdown automata in circular mode. Khrennikov and Yurova
29  modeled the behavior of protein structures using classical automata theory
and investigated the resemblance between the quantum systems and modeling
behavior of proteins. Bhatia and Kumar 30  modeled ribonucleic acid (RNA)
secondary structures using two-way quantum finite automata, which are halted
in linear time. Recently, Duenas-Diez and Perez-Mercader designed molecular
machines for chemical reactions. The native chemical computation has been
implemented beyond the scope of logic gates i.e. with chemical automata 12 .
It has been demonstrated that chemical reactions transcribed in formal
languages, can be recognized by Turing machine without using biochemistry 11 .

##  II Preliminaries

In this section, some preliminaries are given. We assume that the reader is
familiar with the classical automata theory and the concept of quantum
computing; otherwise, reader can refer to the theory of automata 3 , quantum
information and computation 2  7 . Linear algebra is inherited from quantum
mechanics to describe the field of quantum computation. It is a crucial
mathematical tool. It allows us to represent the quantum operations and
quantum states by matrices and vectors, respectively, that obey the rules of
linear algebra. The following are the notions of linear algebra used in
quantum computational theory:

  * •

Vector space (V) 7 : A vector space (V) is defined over the field $\mathbb{F}$
of complex numbers $\mathbb{C}$ consisting of a non-empty set of vectors,
satisfying the following operations:

    * –

Addition: If two vectors $\ket{a}$ and $\ket{b}$ belong to V, then
$\ket{a}+\ket{b}\in V$.

    * –

Multiplication by a scalar: If $\ket{a}$ belongs to V, then $\lambda\ket{a}\in
V$, where $\lambda\in\mathbb{C}$.

  * •

Dirac notation 2 : In quantum mechanics, the Dirac notation is one of the most
peculiarities of linear algebra. The combination of vertical and angle bars
($\ket{}\bra{}$) is used to unfold quantum states. It provides an inner
product of any two vectors. The bra $\bra{b}$ and ket $\ket{a}$ represent the
row vector and column vector, respectively.

$\ket{a}=\begin{bmatrix}\alpha_{1}\\\ \alpha_{2}\\\ \alpha_{3}\\\
\end{bmatrix},\bra{b}=\begin{bmatrix}\beta_{1}^{*}&\beta_{2}^{*}&\beta_{3}^{*}\end{bmatrix},\ket{a}\bra{b}=\begin{bmatrix}\alpha_{1}\beta_{1}^{*}&\alpha_{1}\beta_{2}^{*}&\alpha_{1}\beta_{3}^{*}\\\
\alpha_{2}\beta_{1}^{*}&\alpha_{2}\beta_{2}^{*}&\alpha_{2}\beta_{3}^{*}\\\
\alpha_{3}\beta_{1}^{*}&\alpha_{3}\beta_{2}^{*}&\alpha_{3}\beta_{3}^{*}\\\
\end{bmatrix}$ (1)

where $\beta_{i}^{*}$ indicates the complex conjugate of complex number
$\alpha_{i}$.

  * •

Quantum bit 10 : A quantum bit (qubit) is a unit vector defined over complex
vector space $\mathbb{C}^{2}$. In general, it is represented as a
superposition of two basis states labeled $\ket{0}$ and $\ket{1}$.

$\ket{\phi}=\alpha\ket{0}+\beta\ket{1}$ (2)

The probability of state occurrence $\ket{0}$ is $|\alpha|^{2}$ and $\ket{1}$
is $|\beta|^{2}$. It satisfies that $|\alpha|^{2}+|\beta|^{2}=1$. The two
complex amplitudes ($\alpha$ and $\beta$) are represented by one qubit. Thus,
$2^{n}$ complex amplitudes can be represented by $\it n$ qubits.

  * •

Quantum state 2 : A quantum state $\ket{\psi}$ is defined as a superposition
of classical states,

$\ket{\psi}=\alpha_{1}\ket{w_{1}}+\alpha_{2}\ket{w_{2}}+...+\alpha_{n}\ket{w_{n}}$
(3)

where $\alpha_{i}^{\prime}s$ are complex amplitudes and
$\ket{w_{i}}^{\prime}s$ are classical states for $1\leq i\leq n$. Therefore, a
quantum state $\ket{\psi}$ can be represented as n-dimensional column vector.

$\begin{bmatrix}\alpha_{1}\\\ \alpha_{2}\\\ ...\\\ \alpha_{n}\\\
\end{bmatrix}$ (4)

  * •

Unitary transformation: In quantum mechanics, the transformation between the
quantum systems must be unitary. Consider a state $\ket{\psi}$ of quantum
system at time t:
$\ket{\psi}=\alpha_{1}\ket{w_{1}}+\alpha_{2}\ket{w_{2}}+...+\alpha_{n}\ket{w_{n}}$
transformed into state $\ket{\psi^{\prime}}$ at time t’:
$\ket{\psi^{\prime}}=\alpha_{1}^{\prime}\ket{w_{1}}+\alpha_{2}^{\prime}\ket{w_{2}}+...+\alpha_{n}^{\prime}\ket{w_{n}}$,
where complex amplitudes are associated by
$\ket{\psi^{\prime}(t^{\prime})}=U(t^{\prime}-t)\ket{\psi(t)}$, where U
denotes a time reliant unitary operator satisfies that
$(U(t^{\prime}-t))^{*}U(t^{\prime}-t)=1$, and
$\sum_{i=1}^{n}|\alpha_{i}|^{2}=|\alpha^{\prime}_{i}|^{2}=1$ 2 .

  * •

Hilbert space: A physical system is described by a complex vector space called
Hilbert space $\mathcal{H}$ 7 . It allow us to describe the basis of the
quantum system. The direct sum
$\braket{x}{y}:\mathcal{H}\oplus\mathcal{H}\rightarrow\mathbb{C}$ or inner
product $\braket{x}{v}:\mathcal{H}\otimes\mathcal{H}\rightarrow\mathbb{C}$ of
two subspaces, satisfying the following properties for any vectors:

    * –

Linearity:
($\alpha\bra{x}+\beta\bra{y}$)$\ket{z}=\alpha\braket{x}{z}+\beta\braket{y}{z}$

    * –

Symmetric property: $\braket{x}{y}=\braket{y}{x}$

    * –

Positivity: $\braket{x}{x}\geq 0$ and $\braket{x}{x}=0$ iff $x=0$, where
$x\in\mathcal{H}$.

where $x,y,z\in\mathcal{H}$ and $\alpha,\beta\in\mathbb{C}$.

  * •

Quantum finite automata (QFA) 34 : It is defined as a quadruple
($Q,\Sigma,s_{int},U_{\sigma}$), where

    * –

Q is a set of states,

    * –

$\Sigma$ is an input alphabet,

    * –

Hilbert space $\mathcal{H}$ and $s_{init}\in\mathcal{H}$ is an initial vector
such that $|s_{init}|^{2}=1$,

    * –

$\mathcal{H}_{acc}\subset\mathcal{H}$ and $P_{acc}$ is an acceptance
projection operator on $\mathcal{H}_{acc}$,

    * –

$U_{\sigma}$ denotes a unitary transition matrix for each input symbol
($\sigma\in\Sigma$).

The computation procedure of QFA consists of an input string
$w=\sigma_{n}\sigma_{2}...\sigma_{n}$. The automaton works by reading each
input symbol and their respective unitary matrices are applied on the current
state, starting with an initial state. The quantum language accepted by QFA is
represented as a function $f_{QFA}(w)=|s_{init}U_{w}P_{acc}|^{2}$, where
$U_{w}=U_{\sigma_{1}}U_{\sigma_{2}}...U_{\sigma_{n}}$. The tape head is
allowed to move only in the right direction. Finally, the probability of QFA
being in an acceptance state is observed; i.e. indicating whether the input
string is accepted or rejected by QFA. It is also called as a real-time
quantum finite automaton.

Based on the movement of tape head, QFA is classified as one-way QFA, 1.5-way
QFA, and 2QFA. In 1.5-way QFA, the tape head is permitted to move only in the
right direction or can be stationary, but it cannot move towards the left
direction. It has been proved that it can be designed for non-context free
languages, if the input tape is circular 35 . In this paper, we focused on the
2QFA model due to its more computational power than classical counterparts.

##  III Two-way quantum Finite automata

A quantum finite automaton (QFA) is a quantum variant of a classical finite
automaton. In QFA, quantum transitions are applied on reading the input
symbols from the tape 5 . Two-way quantum finite automaton (2QFA) is a quantum
counterpart of a two-way deterministic finite automaton (2DFA). In 2QFA, the
tape head is allowed to move either in the left or right direction or can be
stationary. The illustration of 2DFA is shown in Fig 1. It consists of a read-
only input tape consisting of a sequence of cells, each one storing a input
symbol. The read-only input tape is scanned by an input head, which is allowed
to move in both directions or can be stationary. At each step of a
computation, a finite state control should be in state from a finite set Q.

Figure 1: Schematic diagram of the “hardware” of 2DFA

######  Definition 1.

_5_ A two-way quantum finite automaton is represented as sextuple
$(Q,\Sigma,\delta,q_{0},Q_{acc},$ $Q_{rej})$, where

  * •

$\it Q$ is a finite set of states.

  * •

$\Sigma$ is an input alphabet,

  * •

Transition function $\delta$ is defined by $\delta:Q\times\Gamma\times Q\times
D\rightarrow\mathbb{C}$, where $\mathbb{C}$ is a complex number,
$\Gamma=\Sigma\cup\\{\\#,\$\\}$ and $D=\\{-1,0,+1\\}$ represent the left,
stationary and right direction of tape head.

  * •

$Q=Q_{acc}\cup Q_{rej}\cup Q_{non}$, where $Q_{non},Q_{acc},Q_{rej}$ represent
the set of non-halting, accepting, rejecting states respectively. The
transition function must satisfies the following conditions:

  1. (i)

Local probability and orthogonality condition:

$\sum_{(q^{\prime},d)\in Q\times
D}^{\forall(q_{1},\sigma_{1}),(q_{2},\sigma_{2})\in
Q\times\Gamma}\overline{\delta(q_{1},\sigma,q^{\prime},d)}\delta(q_{2},\sigma,q^{\prime},d)=\left\\{\begin{array}[]{ll}1~{}~{}q_{1}=q_{2}\\\
0~{}~{}q_{1}\neq q_{2}\\\ \end{array}\right\\}$

  2. (ii)

First separability condition:

$\sum_{q^{\prime}\in Q}^{\forall(q_{1},\sigma_{1}),(q_{2},\sigma_{2})\in
Q\times\Gamma}\overline{\delta(q_{1},\sigma_{1},q^{\prime},+1)}\delta(q_{2},\sigma_{2},q^{\prime},0)+\overline{\delta(q_{1},\sigma_{1},q^{\prime},0)}\delta(q_{2},\sigma_{2},q^{\prime},-1)=0$

  3. (iii)

Second separability condition:

$\sum_{q^{\prime}\in Q}^{\forall(q_{1},\sigma_{1}),(q_{2},\sigma_{2})\in
Q\times\Gamma}\overline{\delta(q_{1},\sigma_{1},q^{\prime},+1)}\delta(q_{2},\sigma_{2},q^{\prime},-1)=0$

For each $\sigma\in\Gamma$, a 2QFA is said to be simplified, if there exists a
unitary linear operator $V_{\sigma}$ on the inner product space such that
$L_{2}\\{Q\\}\rightarrow L_{2}\\{Q\\}$. The transition function is represented
as

$\delta(q,\sigma,q^{\prime},d)=\left\\{\begin{array}[]{l}q^{\prime}V_{\sigma}q\\\
0\end{array}\middle|\;\begin{array}[]{@{}l@{}}\text{if}~{}D(q^{\prime})=d\\\
\text{else}\end{array}\right\\}$ (5)

where $q^{\prime}V_{\sigma}q$ is a coefficient of $\ket{q^{\prime}}$ in
$V_{\sigma}\ket{q}$.

Consider an input string $\it w$, written on the input tape enclosed with both
end-markers such as $\\#w\$$. The computation of 2QFA is progress as follows.
The tape head is above the input symbol $\sigma$ and the automaton is in any
state $\it q$. Then, the state of 2QFA is changed to $q^{\prime}$ with an
amplitude $\delta(q,\sigma,q^{\prime},d)$ and moves the tape head one cell
towards right, stationary and in left direction according to
$\in\\{-1,0,+1\\}$. It corresponds to a unitary evolution in the inner-product
space $\mathcal{H}_{n}$.

A computation of a 2QFA is a chain of superpositions $c_{0},c_{1},c_{2},....,$
where $c_{0}$ denotes an initial configuration. For any $c_{i}$, when the
automaton is observed in a superposition state with an amplitude $\alpha_{c}$,
then it has the form $U_{\delta}\ket{c_{i}}\sum_{c\in
C_{n}}\alpha_{c}\ket{c_{i}}$, where $C_{n}$ represents the set of
configurations. The probability associated with a configuration is calculated
by absolute squares of amplitude. Superposition is said to be valid; if the
sum of the squared moduli of their probability amplitudes is unitary. In
quantum theory, the time evolution is specified by unitary transformations.
Each transition function $\delta$ prompts a transformation operator over the
Hilbert space $\mathcal{H}_{n}$ in linear time.

$U_{\delta}^{w}\ket{q,j}=\sum_{(q^{\prime},d)\in Q\times
D}\delta(q,w(j),q^{\prime},d)\ket{q^{\prime},j+d~{}mod\lvert w\rvert}$

for each $(q,j)\in C_{\lvert w\rvert}$, where $q\in Q,j\in Z_{\lvert w\rvert}$
and extended to $\mathcal{H}_{n}$ by linearity 5 ; 31 .

##  IV Modeling of Chemical Reactions

Figure 2: Representation of language recognition by chemical computation. It
is reproduced from 12  under the Creative Commons CCBY license.

Before, we recognize the chemical reactions using two-way quantum finite
automata model, it is important to show how computational chemistry works. Fig
2. shows the illustration of language recognition by the chemical computation
model. It consists of three parts: (i) a mixed container where the computation
process occurs, (ii) an input translator that translates the chemical aliquots
into input symbols and gives them consecutively depending upon the processing
time, (3) a system to monitor the response of an automaton as a chemical
criterion. Finally, the chemical computation produces well-defined chemical
accept/ reject signatures for the input. For instance, if the number of a’s
and b’s are equal in the input, then the chemical computation produces heat
i.e. an input is said to be accepted. Otherwise, if no heat is released at the
end of computation, then the input is said to be rejected by the system. The
following are the construction of two-way quantum finite state machines of
chemical reactions.

######  Theorem 1.

Two-way quantum finite automata can recognize all regular languages.

###### Proof.

The proof has been shown in 5 . ∎

###  IV.1 Chemical reaction-1 consisting regular language

For an illustrative and visual implementation, we can choose a precipitation
reaction in an aqueous medium such as

$KIO_{3}+AgNO_{3}\rightarrow AgIO_{3}(s)+KNO_{3}$ (6)

If, during computation, a white precipitate of silver iodate is observed, then
the input string is said to be accepted; if the solution is clear from
precipitate, the string has been rejected because there was no reaction.
Therefore, we have chosen the recipes of alphabet symbols a for potassium
iodate ($KIO_{3}$) and b for silver nitrate ($AgNO_{3}$), quantitatively. Fig.
3 shows the chemical representation of symbols a and b, the bimolecular
precipitation reaction 11 . If the precipitate AgIO3 is not presented in the
solution, then the computation is said to be rejected. For example, the input
string $w=aaab$ is said to be accepted due to presence of precipitate or,
equally, the heat has been determined during computation. But, the input
$w=aa$ is said to be rejected due to absence of precipitate, or, precisely,
the heat has not been observed. The Kleene star ($\Sigma^{*}$) operator is a
set of infinite strings of all lengths over input alphabet as well as empty
string ($\epsilon$). Fig. 4 shows the corresponding theoretical 2QFA state
transition graph to recognize $L_{1}$.

Figure 3: Illustration of acid/base reaction of $L_{1}$ Figure 4: State
transition diagram of $L_{1}$

######  Theorem 2.

A language $L_{1}=\\{(a+b)^{*}a(a+b)^{*}b(a+b)^{*}aa^{*}bb^{*}\\}$
representing precipitation reaction in Eq. (6) can be recognized by 2QFA.

###### Proof.

The idea of this proof is as follows. The initial state $q_{0}$ reads a right-
end marker $\\#$ and moves the head towards the right direction. If there is
no occurrence of symbol b, then it shows no precipitate, and the input is said
to be rejected by the 2QFA. Similarly, on reading the symbol b, the state
$q_{0}$ is changed into $q_{1}$. If there is no occurrence of symbol a, then
the state is transformed into rejecting state $q_{r_{1}}$. If the input string
$w\in L_{1}$ contains at least one a and one b, then silver iodate is present
during computation and it is said to be recognized by 2QFA. A 2QFA for $L_{1}$
is defined as follows:

$M_{2QFA}=(Q,\Sigma,q_{0},Q_{acc},Q_{rej},\delta),$

where

  * –

$Q=\\{q_{0},q_{1},q_{2},q_{3},q_{4},q_{5},q_{6},q_{7},q_{a_{1}},q_{a_{2}},q_{r_{1}},q_{r_{2}}\\},$
where $q_{0}$ and $q_{2}$ are used to move the head towards the $ on reading
a’s and b’s respectively. The states $q_{1}$ and $q_{3}$ are used to confirm
the last symbol read by head is a and b, respectively. .

  * –

$\Sigma=\\{a,b\\},q_{0}$ is an initial state,
$Q_{acc}=\\{q_{a_{1}},q_{a_{2}}\\}$ and $Q_{rej}=\\{q_{r_{1}},q_{r_{2}}\\}$.

  * –

The specification of transition functions are given in Table 1.

It can be noted that in 2QFA where transition matrices consist 0 and 1, i.e.
basically a two-way reversible finite automata (2RFA). Therefore, 2QFA can be
designed for all the languages accepted by 2RFA. In transition matrix, each
column and row have exactly only one entry 1. Hence, the dot product of any
two rows is equal to zero. It is known that the language recognition power of
2RFA is an equivalent to 2DFA.

Table 1: Details of the transition functions and head function for $L_{1}$ $V_{\\#}\ket{q_{0}}=\ket{q_{0}}$ |  $V_{a}\ket{q_{0}}=\ket{q_{0}}$ |  $V_{b}\ket{q_{0}}=\ket{q_{1}}$ |  $V_{\\#}\ket{q_{1}}=\ket{q_{2}}$
---|---|---|---
$V_{\$}\ket{q_{0}}=\ket{q_{7}}$ |  $V_{a}\ket{q_{1}}=\ket{q_{2}}$ |  $V_{b}\ket{q_{1}}=\ket{q_{2}}$ |  $V_{b}\ket{q_{2}}=\ket{q_{2}}$
$V_{\$}\ket{q_{7}}=\ket{q_{r_{2}}}$ |  $V_{a}\ket{q_{7}}=\ket{q_{7}}$ |  $V_{a}\ket{q_{2}}=\ket{q_{3}}$ |  $V_{a}\ket{q_{3}}=\ket{q_{1}}$
$V_{\$}\ket{q_{2}}=\ket{q_{5}}$ |  $V_{a}\ket{q_{4}}=\ket{q_{4}}$ |  $V_{b}\ket{q_{3}}=\ket{q_{4}}$ |  $V_{b}\ket{q_{4}}=\ket{q_{3}}$
$V_{\\#}\ket{q_{5}}=\ket{q_{r_{1}}}$ |  $V_{a}\ket{q_{5}}=\ket{q_{6}}$ |  $V_{b}\ket{q_{5}}=\ket{q_{5}}$ |  $V_{b}\ket{q_{6}}=\ket{q_{6}}$
$V_{\$}\ket{q_{4}}=\ket{q_{a_{1}}}$ |  $V_{a}\ket{q_{6}}=\ket{q_{a_{2}}}$ |  |
Head Functions:
$D(q_{0})=(+1),D(q_{1})=(-1),D(q_{2})=(+1),D(q_{3})=(-1),$
$D(q_{4})=(+1),D(q_{5})=(-1),D(q_{6})=(+1),D(q_{7})=(-1),$
$D(q_{a_{1}})=D(q_{a_{2}})=(0),D(q_{r_{1}})=D(q_{r_{2}})=(0)$

∎

###  IV.2 Chemical reaction-2 consisting context-free language

Next, we have considered the context-free language from Chomsky hierarchy
satisfying the balanced chemical reaction between NaOH and malonic acid is as
follows:

$H_{2}C_{3}H_{2}O_{4}+2NaOH\rightarrow Na_{2}C_{3}H_{2}O_{4}+2H_{2}O$ (7)

Figure 5: Illustration of the acid/base reaction of $L_{2}$

The language generated by the above-mentioned chemical reaction is $L_{2}$
consisting Dyck language of all words with balanced parentheses. 2QFA is
designed for $L_{2}$ is as follows:

Figure 6: State transition diagram of $L_{2}$

######  Theorem 3.

A language $L_{2}$ consisting Dyck language of all words with balanced
parentheses can be recognized by 2QFA with probability 1, otherwise rejected
with probability at least $1-\dfrac{1}{N}$, where N is any positive number.

###### Proof.

The idea of this proof is as follows. It consists of three phases. First, the
initial state $q_{0}$ reads a first symbol and both heads start moving towards
the right-end marker $. If the input string starts with closed parentheses,
then it is said to be rejected. On reading the left-end marker $\\#$, the
computation is split into N paths, denoted by $q_{1,0},q_{2,0},...,q_{N,0}$.
Each path possesses equal amplitude $\dfrac{1}{\sqrt{N}}$. Along the N
different paths, each path moves deterministically to the right-end marker
$\$$. Each computational path keep track of the open parentheses with respect
to the closed parentheses. If at the end of computation, the excess of open
parentheses is observed, then it is said to be rejected. It means pH value is
greater than midpoint pH and intermediate gray tone is observed. Secondly, if
there is an excess of closed parentheses, then the darkest gray tone is
observed i.e. pH value is less than midpoint pH. It is said to be rejected by
2QFA with probability $1-\dfrac{1}{N}$. If there is a balanced occurrence of
open and closed parentheses, the input string is said to be accepted with
probability 1. Hence, pH value is equal to midpoint pH and the lightest gray
tone is observed at the end of computation. A 2QFA for $L_{2}$ is defined as
follows:

$M_{2QFA}=(Q,\Sigma,q_{0},Q_{acc},Q_{rej},\delta),$

where

  * –

$Q=\\{q_{0},q_{1},q_{2},q_{3}\\}\cup\\{q_{i,j}\mid 1\leq i\leq N,0\leq
j\leq\text{max}(i,N-i+1)\\}\cup\\{p_{k}\mid 1\leq k\leq N\\}\cup$
$\\{s_{i,0},w_{i,0},r_{i,0}\mid 1\leq i\leq N\\}\cup\\{q_{acc},q_{rej}\\}$,
where $q_{1}$ is used to check whether the first symbol is an open parentheses
or not, $q_{2}$ and $q_{3}$ are used to traverse the input string.

  * –

$\Sigma=\\{(,)\\},q_{0}$ is an initial state, $Q_{acc}=\\{p_{N}\\}$ and
$Q_{rej}=\\{q_{r}\\}\cup\\{p_{k}\mid 1\leq k<N\\}\cup$ $\\{r_{i,0}\mid 1\leq
i\leq N\\}$.

  * –

The specification of transition functions is given in Table 2.

∎

Table 2: Details of the transition functions and head function for $L_{2}$
$V_{\\#}\ket{q_{0}}=\ket{q_{0}},V_{(}\ket{q_{0}}=\ket{q_{0}},V_{)}\ket{q_{0}}=\ket{q_{1}},V_{\\#}\ket{q_{1}}=\ket{q_{r}}$
---
$V_{(}\ket{q_{1}}=\ket{q_{2}},V_{)}\ket{q_{2}}=\ket{q_{2}},V_{(}\ket{q_{2}}=\ket{q_{3}},V_{)}\ket{q_{3}}=\ket{q_{0}}$
$V_{\\#}\ket{q_{2}}=\dfrac{1}{\sqrt{N}}\sum_{i=1}^{N}{\ket{q_{i,0}}}$
$V_{(}\ket{q_{i,0}}=\ket{q_{i,i}},~{}~{}V_{)}\ket{q_{i,0}}=\ket{q_{i,N-i+1}},\text{for}~{}1\leq
i\leq N,V_{(}\ket{q_{i,j}}=\ket{q_{i,j-1}},\text{for}~{}1\leq j\leq i$
$V_{)}\ket{q_{i,j}}=\ket{q_{i,j-1}},\text{for}~{}1\leq j\leq N-i+1,1\leq i\leq
N$
$V_{\$}\ket{q_{i,0}}=\ket{s_{i,0}},V_{)}\ket{s_{i,0}}=\ket{w_{i,0}},V_{(}\ket{s_{i,0}}=\ket{r_{i,0}},\text{for}~{}1\leq
i\leq N$
$V_{\$}\ket{w_{i,0}}=\dfrac{1}{\sqrt{N}}\sum_{k=1}^{N}{exp\left(\dfrac{2\pi
i}{N}ki\right)\ket{p_{k}}},\text{for}~{}1\leq i\leq N$
Head Functions:
$D(q_{0})=(+1),D(q_{1})=(-1),D(q_{2})=(+1),D(q_{3})=(-1),D(q_{i,0})=(+1),\text{for}~{}1\leq
i\leq N$
$D(q_{i,j})=(0),\text{for}~{}1\leq i\leq N,j\neq
0,D(q_{r})=(0),D(p_{k})=(0),\text{for}~{}1\leq k\leq N$
$D(s_{i,0})=(-1),D(w_{i,0})=(+1),D(r_{i,0})=(0),\text{for}~{}1\leq i\leq N$

###  IV.3 Chemical reaction-3 consisting context-sensitive language

To implement a chemical 2QFA for context-sensitive language, we have used
Belousov-Zhabotinsky (BZ) reaction network for the non-linear oscillatory
chemistry 11 , which consists temporal oscillation in the sodium bromate and
malonic acid system 33  as

$3BrO_{3}-+5CH_{2}(COOH)_{2}+3H^{+}\rightarrow
3BrCH(COOH)_{2}+4CO_{2}+2HCOOH+5H_{2}O$ (8)

Figure 7: Illustration of the acid/base reaction of $L_{3}$

In 2019, Duenas-Diez and Perez-Mercader 11  designed chemical Turing machine
for BZ reaction network. The chemical reaction is fed sequentially to the
reactor as $\\{(BrO_{3}^{-})^{n}(MA)^{n}(NaOH)^{n}\\}$, where n$>$0\. It is
transcribed in formal language as $L_{3}=\\{a^{n}b^{n}c^{n}|n>0\\}$. The
symbol a is interpreted as a fraction of sodium bromate, b is used for malonic
acid and symbol c is transcribed as a quantity of NaOH. It is known that
$L_{3}$ is a context-sensitive language and cannot be recognized by finite
automata or pushdown automata with a stack. Although, it can be recognized by
two-stack PDA. We have shown that $L_{3}$ can be recognized by 2QFA without
using any external aid.

######  Theorem 4.

A language $L_{3}=\\{a^{n}b^{n}c^{n}|n>0\\}$ can be recognized by 2QFA in
linear time. For a language $L_{3}=\\{a^{n}b^{n}c^{n}|n>0\\}$ and for
arbitrary N-computational paths, there exists a 2QFA such that for $w\in
L_{3}$, it accepts w with bounded error $\epsilon$ and rejects $w\notin L_{3}$
with probability at least $1-\dfrac{1}{N}$.

Figure 8: State transition diagram of $L_{3}$

###### Proof.

The design of proof for BZ reaction network is as follows. It consists of two
phases. First, the 2QFA traverse the input to check the form
$a^{+}b^{+}c^{+}$. On reading the right-end marker $, the computation is split
into N paths such that $w_{1,0},w_{2,0},...,w_{N,0}$. Secondly, the first path
is used to check the number of $b$’s and $c$’s are equal or not. The second
path is used to check the initial part of an input string to identify it is in
$\\{a^{n}b^{n}\mid n>0\\}$. On reading the right-end marker $, the both paths
are split into N different paths with an equal amplitude
$\dfrac{1}{\sqrt{N}}$. Finally, upon reading the right-end marker #, if the
number of a’s and b’s and b’s and c’s are equal in respective computational
paths. Then, all paths come into N-way Quantum Fourier transform (QFT) and
either one acceptance state or rejectance states are observed. Suppose, if the
input string is not in the corrected form, then all computation paths read the
# at different times. Thus, their amplitudes do not cancel each other and the
input string is said to be rejected with probability $1-\dfrac{1}{N}$.
Otherwise, the input string is said to be recognized by 2QFA with probability
1. ∎

##  V Summary

In summary, 2QFA model can be efficiently designed for balanced chemical
reaction and Belousov-Zhabotinsky (BZ) reaction network with one-sided error
bound, which are halted in linear time. Table 3 shows the language recognition
ability of different computational models. The classical 2DFA and 2PFA are
known to be equal in computational power to one-way deterministic finite
automata (1DFA) 45 ; 46 . It has been proved that 2PFA’s can be designed for
non-regular languages in expected polynomial time. Additionally, it has been
demonstrated that the chemical PDA can be designed for aforementioned chemical
reactions with multiple stacks. The recognition of languages by native
chemical automata cab be found in 12 ; 13 ; 14 . But, we have shown that 2QFA
can recognize such chemical reactions, without any external aid. It has been
proved that 2QFA is more powerful than classical variants, because it follows
the quantum superposition principle to be in more than one state at a time on
the input tape. For the execution, it needs atleast $O(log~{}n)$ quantum
states to store the position of tape head, where n denotes the length of an
input string.

Table 3: Comparison of computational power of models

Languages | Class | 2DFA/2PFA | Chemical FA/PDA | 2QFA
---|---|---|---|---
$L_{1}=\\{(a+b)^{*}a(a+b)^{*}b$ $(a+b)^{*}aa^{*}bb^{*}\\}$ | RL | ✓ | ✓ | ✓
A language $L_{2}$ consisting Dyck language of balanced parentheses | CFL | ✗ |  ✓ (with 1-stack PDA) | ✓
$L_{3}=\\{a^{n}b^{n}c^{n}|n>0\\}$ | CSL | ✗ |  ✓ (with 2-stack PDA) | ✓

  * •

RL, CFL and CSL stand for regular languages, context-free languages and
context-sensitive languages, respectively.

##  VI Conclusion

The enhancement in many existing computational approaches provides momentum to
molecular and quantum simulations at the atomic level. It helps to test new
abstract approaches for considering molecules and matter. Previous attempts to
model the aforementioned chemical reactions used finite automata and pushdown
automata with multiple stacks. In this paper, we focused on well-known
languages of Chomsky hierarchy and modeled them using two-way quantum finite
automata. The crucial advantage of the quantum approach is that these chemical
reactions transcribed in formal languages can be parsed in linear time,
without using any external aid. We have shown that two-way quantum automata
are more superior than its classical variants by using quantum transitions. To
the best of our knowledge, no such modeling of chemical reactions is performed
using quantum automata theory so far. For the future purpose, we will try to
represent complex chemical reactions in formal languages and model them using
other quantum computational models.

##  VII Additional Information

Conflict of interest The authors declare that they have no conflict of
interest.

## Acknowledgement

S.Z. acknowledges support in part from the National Natural Science Foundation
of China (Nos. 61602532), the Natural Science Foundation of Guangdong Province
of China (No. 2017A030313378), and the Science and Technology Program of
Guangzhou City of China (No. 201707010194).

## References

  * (1) M. Conrad, Information processing in molecular systems, BioSystems 5 (1) (1972) 1–14.
  * (2) Y. Benenson, Biocomputers: from test tubes to live cells, Molecular BioSystems 5 (7) (2009) 675–685.
  * (3) S. J. Prohaska, P. F. Stadler, D. C. Krakauer, Innovation in gene regulation: The case of chromatin computation, Journal of theoretical biology 265 (1) (2010) 27–44.
  * (4) B. Bryant, Chromatin computation, PloS one 7 (5) (2012) e35703.
  * (5) M. A. Nielsen, I. Chuang, Quantum computation and quantum information., AAPT, 2002\.
  * (6) J. E. Hopcroft, R. Motwani, J. D. Ullman, Automata theory, languages, and computation, International Edition 24 (2006) 19.
  * (7) P. W. Shor, Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer, SIAM review 41 (2) (1999) 303–332.
  * (8) A. Kondacs, J. Watrous, On the power of quantum finite state automata, in: Proceedings 38th Annual Symposium on Foundations of Computer Science, IEEE, 1997, pp. 66–75.
  * (9) C. Moore, J. P. Crutchfield, Quantum automata and quantum grammars, Theoretical Computer Science 237 (1-2) (2000) 275–306.
  * (10) A. S. Bhatia, A. Kumar, Quantifying matrix product state, Quantum Information Processing 17 (3) (2018) 41.
  * (11) L. Li, D. Qiu, Determination of equivalence between quantum sequential machines, Theoretical Computer Science 358 (1) (2006) 65–74.
  * (12) A. S. Bhatia, A. Kumar, On the power of quantum queue automata in real-time, arXiv preprint:1810.12095 (2018).
  * (13) A. Ambainis, J. Watrous, Two-way finite automata with quantum and classical states, Theoretical Computer Science 287 (1) (2002) 299–311.
  * (14) S. Zheng, D. Qiu, L. Li, J. Gruska, One-way finite automata with quantum and classical states, in: Languages Alive, Springer, 2012, pp. 273–290.
  * (15) A. S. Bhatia, A. Kumar, Quantum $\omega$-automata over infinite words and their relationships, International Journal of Theoretical Physics 58 (3) (2019) 878–889.
  * (16) C. Mereghetti, B. Palano, G. Pighizzini, Note on the succinctness of deterministic, nondeterministic, probabilistic and quantum finite automata, RAIRO-Theoretical Informatics and Applications 35 (5) (2001) 477–490.
  * (17) S. Zheng, D. Qiu, J. Gruska, L. Li, P. Mateus, et al., State succinctness of two-way finite automata with quantum and classical states., Theor. Comput. Sci. 499 (2013) 98–112.
  * (18) A. Yakaryilmaz, A. Say, Succinctness of two-way probabilistic and quantum finite automata, arXiv preprint arXiv:0903.0050 (2009).
  * (19) H. Nishimura, T. Yamakami, An application of quantum finite automata to interactive proof systems, Journal of Computer and System Sciences 75 (4) (2009) 255–269.
  * (20) H. Nishimura, T. Yamakami, Interactive proofs with quantum finite automata, Theoretical Computer Science 568 (2015) 1–18.
  * (21) J. C. Shepherdson, The reduction of two-way automata to one-way automata, IBM Journal of Research and Development 3 (2) (1959) 198–200.
  * (22) S. Zheng, J. Gruska, D. Qiu, On the state complexity of semi-quantum finite automata, RAIRO-Theoretical Informatics and Applications 48 (2) (2014) 187–207.
  * (23) A. Bertoni, C. Mereghetti, B. Palano, Quantum computing: 1-way quantum automata, in: International conference on developments in language theory, Springer, 2003, pp. 1–20.
  * (24) A. S. Bhatia, A. Kumar, On relation between linear temporal logic and quantum finite automata, Journal of Logic, Language and Information (2019) 1–12.
  * (25) D. Qiu, L. Li, An overview of quantum computation models: quantum automata, Frontiers of Computer Science in China 2 (2) (2008) 193–207.
  * (26) A. Ambainis, A. Yakaryılmaz, Automata and quantum computing, arXiv preprint:1507.01988 (2015).
  * (27) A. S. Bhatia, A. Kumar, Quantum finite automata: survey, status and research directions, arXiv preprint:1901.07992 (2019).
  * (28) M. Duenas-Diez, J. Perez-Mercader, Native chemical automata and the thermodynamic interpretation of their experimental accept/reject responses, arXiv preprint:1903.03827 (2019).
  * (29) F. Okubo, T. Yokomori, The computing power of determinism and reversibility in chemical reaction automata, in: Reversibility and Universality, Springer, 2018, pp. 279–298.
  * (30) M. Dueñas-Díez, J. Pérez-Mercader, How Chemistry Computes: Language Recognition by Non-Biochemical Chemical Automata. From Finite Automata to Turing Machines, iScience 19 (2019) 514–526.
  * (31) R. Freivalds, Probabilistic two-way machines, in: International Symposium on Mathematical Foundations of Computer Science, Springer, 1981, pp. 33–45.
  * (32) C. Dwork, L. Stockmeyer, A time complexity gap for two-way probabilistic finite-state automata, SIAM Journal on Computing 19 (6) (1990) 1011–1023.
  * (33) M. Okamoto, T. Sakai, K. Hayashi, Switching mechanism of a cyclic enzyme system: Role as a “chemical diode”, Biosystems 21 (1) (1987) 1–11.
  * (34) A. Hjelmfelt, E. D. Weinberger, J. Ross, Chemical implementation of neural networks and Turing machines, Proceedings of the National Academy of Sciences 88 (24) (1991) 10983–10987.
  * (35) A. Tóth, K. Showalter, Logic gates in excitable media, The Journal of chemical physics 103 (6) (1995) 2058–2066.
  * (36) M. O. Magnasco, Chemical kinetics is Turing universal, Physical Review Letters 78 (6) (1997) 1190.
  * (37) A. Adamatzky, B. D. L. Costello, Experimental logical gates in a reaction-diffusion medium: The xor gate and beyond, Physical Review E 66 (4) (2002) 046112.
  * (38) J. Gorecki, K. Yoshikawa, Y. Igarashi, On chemical reactors that can count, The Journal of Physical Chemistry A 107 (10) (2003) 1664–1669.
  * (39) L. M. Adleman, Molecular computation of solutions to combinatorial problems, Science (1994) 1021–1024.
  * (40) L. Kuppusamy, A. Mahendran, Modelling DNA and RNA secondary structures using matrix insertion–deletion systems, International Journal of Applied Mathematics and Computer Science 26 (1) (2016) 245–258.
  * (41) L. Kuppusamy, A. Mahendran, S. N. Krishna, Matrix insertion-deletion systems for bio-molecular structures, in: International Conference on Distributed Computing and Internet Technology, Springer, 2011, pp. 301–312.
  * (42) T. Krasiński, S. Sakowski, T. Poplawski, Towards an autonomous multistate biomolecular devices built on DNA, in: 2014 Sixth World Congress on Nature and Biologically Inspired Computing (NaBIC 2014), IEEE, 2014, pp. 23–28.
  * (43) A. Khrennikov, E. Yurova, Automaton model of protein: dynamics of conformational and functional states, Progress in biophysics and molecular biology 130 (2017) 2–14.
  * (44) A. S. Bhatia, A. Kumar, Modeling of RNA secondary structures using two-way quantum finite automata, Chaos, Solitons & Fractals 116 (2018) 332–339.
  * (45) J. Wang, Handbook of Finite State Based Models and Applications, CRC press, 2012\.
  * (46) I. Dzelme-Bērziņa, Quantum finite automata and logic, Ph.D. thesis, University of Latvia, Riga (2010).
  * (47) M. Amano, K. Iwama, Undecidability on quantum finite automata, in: STOC ’99: Proceedings of the 31st Annual ACM Symposium on Theory of Computing, ACM, Georgia, USA, 1999, pp. 368–375.
  * (48) A. S. Bhatia, A. Kumar, On the power of two-way multihead quantum finite automata, RAIRO-Theoretical Informatics and Applications 53 (1-2) (2019) 19–35.
  * (49) R. J. Field, E. Koros, R. M. Noyes, Oscillations in chemical systems. ii. thorough analysis of temporal oscillation in the bromate-cerium-malonic acid system, Journal of the American Chemical Society 94 (25) (1972) 8649–8664.
  * (50) M. O. Rabin, D. Scott, Finite automata and their decision problems, IBM journal of research and development 3 (2) (1959) 114–125.