The structure of the generalized Vaidya spacetime containing the eternal naked
singularity

Vitalii Vertogradov

Physics department, Herzen state Pedagogical University of Russia,

48 Moika Emb., Saint Petersburg 191186, Russia

SPb branch of SAO RAS, 65 Pulkovskoe Rd, Saint Petersburg 196140, Russia

<EMAIL_ADDRESS>
## Abstract

In this paper the structure of the generalized Vaidya spacetime when the type-
II of the matter field satisfies the equation of the state $P=\rho$ is
investigated. Satisfying all energy conditions we show that this spacetime
contains the ’eternal’ naked singularity. It means that once the singularity
is formed it will never be covered with the apparent horizon. However, in the
case of the apparent horizon formation the resulting object is a white hole.
We also prove that this spacetime contains only null naked singularity.

Generalized Vaidya spacetime; Naked singularity; White Hole; Vectons

## Introduction

At the end of its lifetime cycle the massive star undergoes the continuous
gravitational collapse. Oppenheimer and Snider were the first who were able to
construct the model of gravitational collapse of the dust cloud [1]. The
result of this collapse is a black hole and it was believed that the result of
such a process is inevitably a black hole. There is a singularity at the
centre of the spherically symmetric black hole. According to the cosmic
censorship conjecture [CCC] any singularity must be hidden under the event
horizon. However, Papapetroo [2] was one of the first who showed that Vaidya
spacetime [3]contains the naked singularity. So Vaidya spacetime is one of the
earliest counterexamples of CCC violation. After that it was understood that
the naked singularity might be the result of the continuous gravitational
collapse [4] (For more detail information see [5]).

When we speak about the result of gravitational collapse, then we should
realize that it might be not only a black hole but also the naked singularity.
The result of such a process depends upon the initial data. One should
understand that this nakedness is temporary and within the short period of
time the singularity will be covered by the apparent horizon. Hence there is a
question about the existence of the eternal naked singularities. Under the
notion ’eternal’ we assume that the singularity was formed, but will never be
covered with the apparent horizon. It was shown that in the case of the
gravitational collapse of the scalar field [6] and the generalized Vaidya
spacetime [10] the result might be the eternal naked singularity formation.
However, in these models one has to violate the strong energy condition [8]
and assumes that the matter possesses the negative pressure. So there is a
question if we can construct the model of the eternal naked singularity
formation without violating any energy condition.

In this paper, we consider the generalized Vaidya spacetime [9] when the type
II of the matter field satisfies the equation of the state $P=\rho$, where $P$
and $\rho$ are pressure and energy density respectively. The type-I of the
matter field represents the null dust. The equation of the state $P=\rho$ is
so-called the extremely rigid equation of the state and it corresponds to the
particle - vecton [12]. We can also notice that the extremely rigid equation
of the state corresponds to the case when the speed of the sound is equal to
the speed of light and often the matter is called as the stiff fluid. In this
article we will show that this space-time contains the eternal naked
singularity and also a white hole. The paper is organized as follows: in sec.
II we consider the generalized Vaidya spacetime when the II-type of the matter
spacetime possesses the extremely rigid equation of the state. In sec. III we
consider the eternal naked singularity formation. In sec. IV white hole
formation is investigated. Sec. V is the conclusion. Throughout the paper
$c=G=1$ system of units will be used.

## Generalized Vaidya spacetime
with the extrimely rigid equation of the state

The generalized Vaidya spacetime has the following form [9]:

$\begin{split}ds^{2}=-\left(1-\frac{2M(r,v)}{r}\right)dv^{2}+2\varepsilon
dvdr+r^{2}d\Omega^{2}\,,\\\ d\Omega^{2}=d\theta^{2}+\sin^{2}\theta
d\varphi^{2}\,,\end{split}$ (3.1)

here $M(v,r)$ \- the mass function depending on coordinates $r$ and $v$ which
corresponds to advanced/retarded time, $\varepsilon=\pm 1$ depending on
ingoing/outgoing shells. We can write down the energy momentum tensor in the
following form:

$T_{ik}=T^{(n)}_{ik}+T^{(m)}_{ik}\,,$ (3.2)

where the term $T^{(m)}_{ik}$ corresponds to the matter field of the type-I
and the other $T_{ik}^{(n)}$ corresponds to the matter field of the type-II
[7]. Now let us write down the expression of the energy momentum tensor: [9]

$\begin{split}T^{(m)}_{ik}=\mu L_{i}L_{k}\,,\\\
T^{(n)}_{ik}=(\rho+P)(L_{i}N_{k}+L_{k}N_{i})+Pg_{ik}\,,\\\
\mu=\frac{2\varepsilon\dot{M}}{r^{2}}\,,\\\
\rho=\frac{2M^{\prime}}{r^{2}}\,,\\\ P=-\frac{M^{\prime\prime}}{r}\,,\\\
L_{i}=\delta^{0}_{i}\,,\\\
N_{i}=\frac{1}{2}\left(1-\frac{2M}{r}\right)\delta^{0}_{i}-\varepsilon\delta^{1}_{i}\,,\\\
L_{i}L^{i}=N_{i}N^{i}=0\,,\\\ L_{i}N^{i}=-1\,.\end{split}$ (3.3)

here $P$ \- pressure, $\rho$ \- density, $\mu$ \- the energy density of the
null dust. And $L,N$ \- two null vectors. This model must be physically
reasonable, so the energy momentum tensor should satisfy weak, strong and
dominant energy conditions. It means that $\rho$ must be positive and for any
non-spacelike vector $v^{i}$:

$T_{ik}v^{i}v^{k}>0\,,$ (3.4)

and the vector $T_{ik}v^{i}$ must be timelike. Strong and weak energy
conditions demand:

$\begin{split}\mu\geq 0\,,\\\ \rho\geq 0\,,\\\ P\geq 0\,.\end{split}$ (3.5)

The dominant energy condition imposes following conditions on the energy
momentum tensor:

$\begin{split}\mu\geq 0\,,\\\ \rho\geq P\geq 0\,.\end{split}$ (3.6)

The properties of the generalized Vaidya spacetime have been studied for the
equation of the state $P=\alpha\rho$ where $\alpha$ belongs to the interval
$(-1,\frac{1}{3}$ in articles [13, 11, 14]. If we satisfy this equation of the
state, then the mass function $M(r,v)$ has the form

$M(r,v)=C(v)+D(v)r^{1-2\alpha}\,,$ (3.7)

where $C(v)$ and $D(v)$ are positive functions of time $v$. One should note
that the speed of sound in the medium depends upon $\alpha$ in the following
way :

$v_{sound}=\sqrt{\alpha}\,.$ (3.8)

We want to study the generalized Vaidya spacetime with the extremely rigid
equation of the state. For this purpose, we must consider the equation of the
state:

$P=\rho\,.$ (3.9)

In this case the mass function has the form (3.7) with $\alpha=1$. One can see
from (3.8) that the speed of sound in this case is equal to the speed of
light. However, we must consider function $D(v)$ properly. Now it can’t be
positive because the density expression is given by:

$\rho=-\frac{D(v)}{r^{4}}<0\,,$ (3.10)

and if $D(v)>0$ then we violate energy conditions (3.5) and (3.6). So to
satisfy these conditions we must assume that $D(v)<0$ or we can write:

$D(v)=-D^{\prime}(v)\,,$ (3.11)

where now function $D^{\prime}(v)$ is positive. We want to investigate the
question about the singularities in this spacetime which are formed at $v=0$
and $r=0$. But in this case we must impose some extra conditions on functions
$C(V)$ and $D^{\prime}(v)$ not to violate energy conditions (3.5) and (3.6)
because if $\dot{C}(v)<\frac{\dot{D^{\prime}}(v)}{r}$ at some point then $\mu$
becomes negative and we violate our energy conditions. So we must consider two
cases:

  1. 1.

$D^{\prime}(v)\equiv\lambda$, where $\lambda$ is positive real constant. In
this case $\dot{D^{\prime}}(v)\equiv 0$ and if we impose the following
condition $C(V)\geq 0$ then we satisfy all energy condition.

  2. 2.

$D^{\prime}(0)=0$ but in this case even if $\lim\limits_{v\to 0,r\to
0}\frac{D^{\prime}(v)}{r}=X_{1}$ where $X_{1}$ finite positive constant and if
$C(v)>X_{1}$ then anyway we would violate energy conditions at the later
stage.

Also $C(v)$ must satisfy the following condition:

$\begin{split}\dot{C}(v)\geq 0\,\mathrm{if}\,\varepsilon=+1\,,\\\
\dot{C}(v)\leq 0\,\mathrm{if}\,\varepsilon=-1\,.\end{split}$ (3.12)

## The structure of the singularity

The apparent horizon equation in the case of generalized Vaidya spacetime with
extremely rigid equation of the state is given by:

$\frac{2C(v)}{r}-\frac{2\mu}{r^{2}}-1=0\,.$ (4.13)

The time of the singularity formation is $v=0$. So solving the above equation,
we obtain:

$\begin{split}R^{2}-2C(v)r+2\mu=0\,,\\\
\frac{D}{2}=C^{2}(v)-2\mu\,.\end{split}$ (4.14)

From this equation we can conclude that if

$C^{2}(v)<2\mu\,,$ (4.15)

then the apparent horizon will be never formed. Now let’s consider the
question about the existence of null radial geodesics which terminate at the
central singularity in the past. The radial null geodesic equation has the
following form:

$\frac{dv}{dr}=\frac{2\varepsilon r^{2}}{r^{2}-2C(v)r+2\mu}\,.$ (4.16)

If we consider the limit $\lim\limits_{v\to 0,r\to 0}\frac{dv}{dr}$ then we
can see that In the denominator of this fraction (4.16) we have a finite
positive number (see (4.14).). And as the result we have:

$\begin{split}X_{0}=\lim\limits_{v\to 0,r\to 0}\frac{dv}{dr}=0\,,\\\
v=\text{const}\,.\end{split}$ (4.17)

Usual and generalized Vaidya spacetimes we obtain from Schwarzschild spacetime
assuming that the time axe is the radial null geodesic in Schwarzschild
spacetime because of it, we can have $X_{0}=0$ as the condition for the
existence of the radial null geodesic. It is also because to obtain the Vaidya
metric one should do the following coordinate transformation of the
Schwarzschild spacetime:

$v=t+\varepsilon(r+2m\ln\left|r-2m\right|)\,.$ (4.18)

So in the case $v=const.$ we have the eternal naked singularity
formation111Under the notion ’eternal’ we understand that the singularity is
formed and will be never covered with the apparent horizon. However, not not
all $v=const.$ suits us, the result depends on $\varepsilon$. If
$\varepsilon=+1$ \- falling matter, then the null adial geodesic $v=const.$ is
past oriented. So to have future oriented radial null geodesic $v=const.$
which originate at the central singularity in the past one should consider
outgoing Vaidya spacetime $\varepsilon=-1$ i.e. the radiating Vaidya
spacetime.

###  § 4.1 The radial timelike geodesics

We have found out that there is a family of the radial null future-directed
geodesics which terminate at the central singularity in the past. Now we are
interested if not only massless particles can move away from the central
singularity but also particles which possess the mass. For this purpose, we
must consider the timelike geodesic equation. We can notice that the
generalized Vaidya spacetime extremely rigid equation of the state doesn’t
depend on the angle $\varphi$. Hence we can conclude that the angular momentum
$L$ is conserved. However, the metric depends on the time $v$ so the energy is
not constant, but some function of the time $v$. To derive the timelike
geodesic equation, we should consider the lagrangian. In this case it has the
following form:

$2\mathcal{L}=-\left(1-\frac{2C(v)-\frac{2\mu}{r}}{r}\right)\dot{v}^{2}+2\dot{v}\dot{r}+r^{2}(\dot{\theta}^{2}+\sin^{2}\theta\dot{\varphi})\,,$
(4.19)

where sign dot denotes the partial derivative on affine parameter $\tau$. We
are interested in the radial geodesics so because of it, we can demand
$\dot{\theta}=\dot{\varphi}=0$. Taking into account this one and (4.19) we can
find the energy $E$ as a function of the time $v$:

$-E(v)=\frac{\partial\mathcal{L}}{\partial\dot{v}}=-\left(1-\frac{2C(v)-\frac{2\mu}{r}}{r}\right)\dot{v}+\dot{r}\,.$
(4.20)

We impose the condition $g_{ik}\dot{x}^{i}\dot{x}^{k}=-1$ to obtain the
timelike geodesic. In our case this condition gives:

$-1=(-E+\dot{r})\dot{v}\,.$ (4.21)

Now if we use (4.20) the equation (4.21) gives:

$\dot{r}=\pm\sqrt{E^{2}(v)-\left(1-\frac{2C(v)-\frac{2\mu}{r}}{r}\right)}\,.$
(4.22)

Now substituting (4.22) into (4.20) we obtain:

$\dot{v}=\frac{-E-\sqrt{E^{2}(v)-\left(1-\frac{2C(v)-\frac{2\mu}{r}}{r}\right)}}{-\left(1-\frac{2C(v)-\frac{2\mu}{r}}{r}\right)}\,,$
(4.23)

we picked up sign ”+” because we are interested in outgoing geodesics. Now
combining (4.22) and (4.23) we finally obtain the timelike geodesic equation:

$\frac{dv}{dr}=\frac{-E-\sqrt{E^{2}(v)-\left(1-\frac{2C(v)-\frac{2\mu}{r}}{r}\right)}}{-\left(1-\frac{2C(v)-\frac{2\mu}{r}}{r}\right)\sqrt{E^{2}(v)-\left(1-\frac{2C(v)-\frac{2\mu}{r}}{r}\right)}}\,.$
(4.24)

Now if we consider the limit $\lim\limits_{v\to 0,r\to 0}\frac{dv}{dr}$ then
we can easily see that this limit is equal to zero. Hence we can conclude that
the equation (4.24) has the solution $v=\text{const.}$. But the case
$v=\text{const.}$ is suited only for the radial null geodesics (not for
timelike) so we can conclude that there exists only the radial null geodesic
which terminates at the central singularity int the past.

## The white hole formation

We can see that (4.14) can have positive solutions if

$C^{2}(v)\geq 2\mu\,.$ (5.25)

So we have solutions:

$r_{\pm}=C(v)\pm\sqrt{C^{2}(v)-2\mu}\,.$ (5.26)

It is worth notecing that these roots are both positive. Let’s consider the
region $0\leq r\leq r_{-}$. The existence of a famile of the radial null
geodesics which terminate at the central singularity in the past we proved in
the previous section. Also, we can notice that the limit $\lim\limits_{r\to
0,v\to 0}\frac{dv}{dr}=0$ i.e. this limit is not negative and it means that we
have the movement only from the central singularity and we don’t have the
movement towards it. Now we should prove that the expansion $\theta$ [8] is
positive everywhere in this region. The expansion in the case of generalized
Vaidya spacetime extremely rigid equation of the state has the following form
[10]:

$\theta=e^{-\gamma}\frac{2}{r^{3}}\left(r^{2}-2C(v)r+2\mu\right)\,,$ (5.27)

where the parameter $\gamma$ doesn’t have any impact on the sign $\theta$. So
only the expression in the parentheses has the impact of the sign $\theta$.
But from (4.14) we can notice that this expression is positive in the region
$0\leq r\leq r_{-}$. We have proved that the expansion $\theta$ is positive in
the region $0\leq r\leq r_{-}$. Hence this region corresponds the white hole
solution.

## Conclussion

In this article we considered the particular solution of the generalized
Vaidya spacetime when the type-II of the matter fields satisfies the extremely
rigid equation of the state$P=\rho$. This solution corresponds to vecton
particles. In this spacetime the speed of sound is equal to the speed of
light. In this case, satisfying all energy conditions, we prove that this
spacetime contains the ’eternal’ null naked singularity. However, this
nakedness is possible only for radiating metric, i.e. $\varepsilon=-1$. The
other possibility is the white hole formation. In comparison with naked
singularity the white hole also contains the family of non-spacelike future-
directed geodesics which terminates at the central singularity in the past.
But as compared to the naked singularity, the white hole singularity is
covered with apparent horizon.

One should note that the case of the stiff fluid is similar to the usual
charged Vaidya solution 222Under ’usual’ we understand that the mass function
depends only on time $v$ but the black hole possesses the charge $Q$.. This
metric doesn’t violate the weak energy condition [15] but one can maintain
this condition due to charge and the Lorentz force. However, here the function
$D(v)\neq Q(v)$ because in the case of the absence of the matter field type-II
$D(v)\equiv 0$. Thus, there is no any Lorentz force which can prevent the
violation of the weak energy condition in the general case because of it we
have to demand the restrictions (3.11).

We constructed the model of the ’eternal’ naked singularity without violating
any energy conditions. However, this singularity is not timelike. To obtain
the timelike ’eternal’ naked singularity we either should consider another the
equation of the state or another spacetime. Obtaining of timilike eternal
naked singularity is the question of future investigations.

acknowledgments The author says thanks to professor Pankaj Joshi for
scientific discussion and grant Num. 22-22-00112 RSF for financial support.

## References

  * [1] J.R. Oppenheimer, H.Snyder, Phys. Rev. 56, 455-459 (1939).
  * [2] A. Papapetrou, In: A Random Walk in Relativity and Cosmology. Wiley Eastern, New Delhi (1985).
  * [3] P. C. Vaidya. Nonstatic solutions of Einstein’s field equations for spheres of fluids radiating energy. Phys. Rev., 83:10, 1951.
  * [4] Pankaj S. Joshi Gravitational collapse and spacetime singularities. Cambridge University Press. 2007.p.273.
  * [5] Pankaj S. Joshi, D. Malafarina, ”Recent development in gravitational collapse and spacetime singularitits”. Arxiv.org/abs/1603.02848.
  * [6] Goswami R. Gravitational collapse of a selfinteracting scalar field / R. Goswami, P.S. Joshi // Modern Physics Letters A Vol. 22, No. 1. 2007. P. 65-74.
  * [7] S. W. Hawking and G. F. R. Ellis, The Large Scale Structure of Space-Time. Cambridge University Press, Cambridge (1973).
  * [8] Poisson E. A Relativist’s Toolkit: The Mathematics of Black-Hole Mechanics. Cambridge University Press (2007).
  * [9] A. Wang, Yu.Wu ”‘Generalized Vaidya solutions”’. Gen Relativ. Gravit., 31, 107 (1999)
  * [10] Vertogradov V. The eternal naked singularity formation in the case of gravitational collapse of generalized Vaidya spacetime Volume No.33, Issue No. 17.
  * [11] V.D. Vertogradov ”‘Naked singularity formation in generalized Vaidya space-time”’. Grav. Cosmol. 22, 220-223 (2016).
  * [12] Ya.B. Zel’dovich, I.D. Novikov Theory of gravitation and the evolution of stars, Moskva: Nauka, 1971, 484 p.
  * [13] Vertogradov V.D. Conditions for the naked singularity formation in generalized Vaidya spacetime. Journal of Physics: Conference Series. 2016. Vol. 769. Art. 012013.
  * [14] Vertogradov V. Gravitational collapse of Vaidya spacetime. International Journal of Modern Physics: Conference Series Vol. 41. 2016. Art. 1660124.
  * [15] Soumyabrata Chatterjee, Suman Ganguli, Amitabh Virmani Charged Vaidya Solution Satisfies Weak Energy Condition General Relativity and Gravitation, 48(7), 1-22, 2016


# Analysis of chaotic dynamical systems with autoencoders

N. Almazova1 G. D. Barmparis1,2,∗<EMAIL_ADDRESS>G. P. Tsironis1,2
1National University of Science and Technology MISiS, Leninsky prosp. 4,
Moscow 119049, Russia  2Institute of Theoretical and Computational Physics and
Department of Physics, University of Crete, P.O. Box 2208, 71003 Heraklion,
Greece

###### Abstract

We focus on chaotic dynamical systems and analyze their time series with the
use of autoencoders, i.e., configurations of neural networks that map
identical output to input. This analysis results in the determination of the
latent space dimension of each system and thus determines the minimal number
of nodes necessary to capture the essential information contained in the
chaotic time series. The constructed chaotic autoencoders generate similar
maximal Lyapunov exponents as the original chaotic systems and thus encompass
their essential dynamical information.

††preprint: AIP/123-QED

> The quantitative study of complex systems may be assisted by the use of
> methods derived from neural networks. The analysis of chaotic time series
> presents an applied challenge towards determining the features involved and
> particularly the effective complexity they include. This work shows that the
> use of chaotic autoencoders made of neural networks leads to essential
> information on the time series and the chaotic systems they arise from. The
> autoencoders we construct map identically the input chaotic time series to
> the output. The process passes through the latent space that acts as an
> information bottleneck for this reconstruction. We estimate the dimension of
> the latent space numerically and find it smaller than possible embedding
> dimensions of the reconstruction. The dimension of the latent space hints
> towards the complexity of the dynamical system underlying the time series.
> Furthermore, we show that the constructed chaotic autoencoders produce
> maximal Lyapunov exponents that are very close to those obtained directly
> from the equations of motion of the chaotic systems. These features show
> that the constructed chaotic autoencoders are faithful representations of
> the time series generating dynamical systems, and through their latent space
> dimension, we obtain an estimate of their complexity.

##  I Introduction

Nonlinear dynamical systems, either continuous or discrete, may produce chaos,
i.e. irregular evolution similar to probabilistic, stochastic motion.
Continuous autonomous chaotic systems are described by at least three coupled
ordinary nonlinear differential equations that, when solved numerically,
result in trajectories with irregular features. In general, non-Hamiltonian
cases, the trajectory may fall on a strange attractor that is a chaotic yet
distinct structure in the system phase space. Understanding the complex
dynamics in phase space is significant especially in the mathematical analysis
of the inverse problem. In the latter we are given a time series of discrete
data and we want to find whether this set derives from a chaotic dynamical
system or it is purely stochastic Anishchenko ; Abarbanel ; Kantz ; Bradley .
This is a significant problem that stems from the need to understand
experimental data as well as other data that are readily available in our
times. In this problem, we are given a time series of data and we are
interested in reconstructing the phase space of the associated dynamical
system Abarbanel2 ; Deyle ; Verdes . The methods that have been developed over
the last thirty or more years work quite well and the area is mature.
Nevertheless, since these methods are necessarily cumbersome and involve
additionally empirical inference it would be interesting to experiment with
the use of methods derived from Artificial Intelligence (AI). More
specifically, Machine Learning (ML) provides with an ensemble of well tested
methods that may assist or even surpass traditional methods of analysis of
chaotic systems Goodfellow ; Jordan . The basic aim of the present work is to
use unsupervised learning methods and attempt to reconstruct the dynamical
features of the data without relying on the traditional methods.

We propose an approach that employees autoencoders, a ML tool that is used
widely, as a dimensionality reduction technique. Autoencoders are implemented
with neural networks and are in effect unit projection operators where the
output is as close as possible to the input. In the process of this unit
projection, data pass through the latent space where they suffer
dimensionality reduction. The effective dimension of the latent space gives
the degree of possible data compression and thus eliminates original dimension
redundancy. Since the neural network uses fewer nodes in latent space than the
number of input nodes we accomplish compression of the original data. This
packing allows detecting the basic information contained in the data. This
compression technique involves data regularization in latent space that is
responsible for the selection of the optimal number of operating nodes. This
strategy allows us to find the essential number of nodes for the trained
network that recapture the complexity of the original time series.

Once the minimal configuration in latent space is found, we need to compare
the chaotic autoencoder with the original system. This is accomplished through
the comparison of the corresponding maximal Lyapunov exponents, i.e. that of
the original trajectories of the system with the one calculated from the
reconstructed trajectories generated through the autoencoder. The Lyapunov
exponents are a general attribute of nonlinear dynamical systems as they allow
to distinguish the local stability and invariant sets in a particular
direction of attractors. The largest Lyapunov exponent detects the mean
divergence of nearest neighbors in two initial nearby trajectories. The
positive value of the exponent defines the degree of chaos in a particular
system. There exist different methods for finding the largest Lyapunov
exponent Wolf ; Sano ; Rosenstein ; Eckmann ; Parlitz . In this work we use
the largest Lyapunov exponent as a measure of the level of chaos resulting
from neural network chaotic autoencoders.

Machine Learning techniques are used in areas far beyond Computer Science.
Neural networks (NN) are applied for a wide range of tasks that permit solving
problems in engineering as well as science while deep learning expands the
ability of generalization, interpolation, etc Jordan ; Nielsen . Artificial
intelligence helps to solve tasks through appropriate data handling but
nevertheless they cannot give deep insight to the solution processes they use.
They operate by allowing the computers to learn the features of data and
subsequently understand the hierarchies involved. Several Machine Learning
techniques, like the reservoir computing Ott ; Budhiraja , the long short-term
memory (LSTM) Vlachas  and the O-LSTM Neofotistos , have been used to
investigate the time evolution of complex systems when the time series is
known Brunton ; Dercole ; Lusch . Autoencoders are also used in the physical
systems due to their main advantage of dimensionality reduction Goodfellow ;
Lopez . An autoencoder is a neural network that maps the input data to the
output Kramer . It presents a method of unsupervised learning with the
compression of data occurring in the latent space that forms a kind of
information bottleneck, due to the smaller number of nodes it includes. The
target data for an autoencoder are the input data. Traditionally, autoencoders
are used as dimensionality reduction techniques and salient features of
learning. The construction of an autoencoder includes two main parts. The
first part is an encoder function, $h(x)$, that compresses the input
information, $x$, while the second part is the decoder, $g(h)$ takes as input
the compressed information and by decompressing it reconstructs the initial
input of the encoder. The output function can be written as
$\hat{x}=g\big{(}h(x)\big{)}$, where the input sequence should be learned
“perfectly". In other words, the autoencoder attempts to copy the data to the
output while first they pass through a compression bottleneck. The learning
results are obtained through a minimization of the appropriate loss function.

In this article we focus on the construction of chaotic autoencoders. In the
next section, we detail features of dynamical systems such as phase space
reconstruction and largest Lyapunov exponents and report on the methods for
their calculation. Subsequently, we define the chaotic autoencoders for three
prototypical dynamical systems and show how the information compression
operates. We then show how the Lyapunov exponents are evaluated in the chaotic
autoencoders, and finally, we conclude.

##  II Description of the chaotic systems

A deterministic chaotic signal behaves as something intermediate between a
regular and a stochastic motion. The nonlinear aspect of this behavior of the
signal makes difficult the prediction and classification of the physical
system, from which the motion comes from Deyle ; Sauer ; Hilborn ; Lorenz . A
chaotic motion produces specific structures in phase space. Systems of
differential equations with three or more degrees of freedom may show chaos in
their evolution in real or phase space. The number of degrees of freedom of a
system describes the number of the necessary independent variables to specify
the dynamical state of the system. Chaotic systems are described through their
time-evolution equations, the values of their parameters, and the selected
initial conditions. We may express the ordinary differential equations guiding
the evolution of a state using a variable
$x(t)=[x_{1}(t),x_{2}(t),\ldots,x_{i}(t)]$ where $x_{i}(t)$ is the variable
describing the $i-th$ degree of freedom at time, $t$, as:

$\frac{dx(t)}{dt}=f(x(t))$ (1)

where, $f$, describes the nonlinear function that characterizes the equations
of motion of the system. Given the system of differential equations guiding a
dynamical system, we can solve them numerically and follow its evolution in
time. Many dynamical systems have been investigated using this methodology.
Here, we focus on prototypical systems well known in the literature of chaotic
dynamics, viz. the Rössler system and the Lorenz system. Since we want to have
a more general idea of the effectiveness of our method, we use both the
Lorenz63 system that contains three degrees of freedom and the Lorenz96 one
with multiple dependent variables and external forcing. We will describe these
systems below.

###  II.1 Delay-coordinate embedding

A physical system can only be measured discretely in time, with a given time
step, and thus it can be represented as a time series of data. In the case of
a chaotic system, we can use several methods to determine the features of the
chaotic data. One of these methods is the time-delay embedding or
reconstruction scheme Abarbanel2 ; Hilborn ; Broomhead . This method allows
for the recovery of the entire multidimensional observation in state-space
based on the time series of just one of its variables. Let us assume that a
time series, $X$, represents a trajectory of a physical system recorded with a
fixed time step as an $N$-points time series. The reconstructed trajectory,
$X$, can be represented as an $(M\times m)$ matrix,

$X=\big{[}X_{1},X_{2},\ldots,X_{M}\big{]}^{T},$ (2)

where, $M=N-(m-1)\tau$, $m$ is an embedding dimension, $\tau$ is a time lag
(delay) and $X_{i}$ is a vector representing the state of the system at
discrete time $i$, given by:

$X_{i}=\big{[}x_{i},x_{i+\tau},\ldots,x_{i+(m-1)\tau}\big{]},$ (3)

with, $x_{i}$, being the, $i-th$, data point of the time series. In order to
find the embedding dimension for a dynamical system, one can implement the
Takens theorem Takens . Takens introduced the idea that a complete part of
dynamics with basic dimension $d_{s}$ might be observed in an embedding space
if the latter has $2d_{s}+1$ or higher dimensions. There are a number of
practical implementations of Takens theorem although these methods are not
necessarily universal Kantz ; Packard ; Rosenstein2 . In the present work, we
use a time lag equal to one and an appropriate embedding dimension that
satisfies the Takens theorem.

##  III Methodology

###  III.1 Data Preparation

In this work, we analyze three chaotic systems, viz. that of Rössler,
Lorenz63, and Lorenz96. For the sake of brevity, we present the details of
each system in Table 1. We integrate the differential equations of each system
by means of the 4th order Runge-Kutta method, using a fixed time step equal to
0.005, for 300000 time steps. The first 50000 time steps of each solution were
removed as transient points. We use the x-coordinate of the Rössler and the
Lorenz63 systems, and the 19th coordinate (out of forty) of the Lorenz96 as
the training sequences. The 250000 time steps of each training sequence are
normalized between zero and one and split into training (80%) and testing
(20%) sets. The training and testing sets are split into smaller pieces using
a sliding window of a given size, $W$.

System | Equations | Parameters
---|---|---
Rössler |  |
| $\frac{dx}{dt}=y-z$ | a = 0.1
$\frac{dy}{dt}=x+ay$ | b = 0.1
$\frac{dz}{dt}=b+z(x-c)$ | c = 14
Lorentz63 |  |
| $\frac{dx}{dt}=\sigma(y-x)$ |  $\sigma$ = 10
$\frac{dy}{dt}=x(\rho-z)-y$ |  $\rho$ = 28
$\frac{dz}{dt}=xy-\beta z$ |  $\beta$ = 8/3
Lorentz96 |  |
|  | N = 40
$\frac{dx_{i}}{dt}=(x_{i+1}-x_{i-2})x_{i-1}-x_{i}+F$ |
|  $F$ = 8.15
Table 1: A graphical representation, the equations and the parameters of each
one of the chaotic systems in consideration.

###  III.2 Machine Learning: Setting up the chaotic autoencoder

A stacked autoencoder composed of two symmetric networks, the encoder, $h(x)$,
and the decoder, $g(h)$, is used to reconstruct each one of the chaotic
systems in this work. The encoder consists of an input layer with 30 nodes,
i.e., equal to the window size, $W$, used to split the input sequences. Then,
two fully connected hidden layers of 22 and 15 nodes with a Sigmoid and ReLU
activation function, respectively, guide the model to the latent space. The
latent space is represented by a fully connected layer of 10 nodes and a ReLU
activation function. L1 regularization (Lasso Regression) is applied to the
output of the nodes of the latent space to force nodes with less information
to zero. The L1 regularization applied to the output of a layer adds to the
loss function a value proportional to the regularization parameter, $\alpha$,
to shrink nodes with less information. The regularization parameter, $\alpha$,
is used as a hyper-parameter of the model.The larger the regularization
parameter, the more we constrain the autoencoder and thus we increase the
error in the reconstructed sequence. The smaller the regularization parameter,
the larger the active part of the latent space and thus we decrease the error
in the reconstructed sequence. After finishing with the training process, we
statistically analyze the contribution of the number of active nodes in the
latent space to accuracy of the the reconstructed sequence. Nodes that output
values smaller than the regularization parameter contain non-important
information and their contribution to the reconstructed sequence is
negligible. Nodes with values larger than the value of $\alpha$ are counted as
non-zero nodes with significant contribution to the output. The decoder
network is constructed as the mirror image of the encoder one. The autoencoder
model is trained for 7500 epochs with a batch size equal to 32, using the Adam
stochastic optimization method (ADAM, ) with a learning rate of 0.001, while
all its weights are initialized using the Xavier uniform initializer Xavier ,
as implemented in Tensorflow/Keras TF ; Keras . In Fig. (1), we present a
schematic representation of the proposed model.

Figure 1: A schematic representation of the proposed model. The trajectories
of a given attractor are transformed into input sequences, $x$, that pass
through an encoder model, $h(x)$, get compressed into the latent space and
pass through the decoder, $g(h)$, to reconstruct the original information,
$\hat{x}=g(h(x))$.

As mentioned in the Introduction section, an autoencoder is a neural network
that acts as an identity operator and maps the input data ideally to identical
output data. In performing this operation through the latent space, it must
learn through minimizing a loss function. The simplest example of a loss
function is that of the mean squared error (MSE) , i.e.

$MSE(x,\hat{x})=\frac{1}{W}\sum_{k=1}^{W}\big{(}x_{k}-\hat{x}_{k}\big{)}^{2},$
(4)

where, $W$ is the length of the input sequence, $x$, represents the real -
input data (a set of $x_{k}$ real samples), $\hat{x}=g(h(x))$ represents the
reconstructed - output data (a set of $\hat{x}_{k}$ reconstructed samples) of
a general ML system, while the loss function "punishes" the difference between
$x$ and $\hat{x}$, i.e., the real and reconstructed data. Here, to train the
autoencoder we additionally apply the sparsity penalty (L1 regularization),
$\Omega(h)=\alpha|h(x)|$, on the output of the encoder, $h(x)$, (after the
activation function), where, $\alpha$, is the regularization parameter. Then
the loss function, $L$, is represented as:

$L\big{(}x,\hat{x}\big{)}=MSE+\Omega(h(x)).$ (5)

This feature permits the autoencoder to limit the number of nodes at the
latent space to a minimum required number, that can be used to reconstruct the
original signal correctly. This type of regularization forces the
participation of nodes in latent space to be zero when their contribution to
the output of the latent space is less than the regularization parameter.
Given the fact that we want to investigate the performance of the model as a
function of the regularization parameter, $\alpha$, in order to demonstrate
the effect of the regularization and find the minimum number of nodes that can
correctly reconstruct the input sequences, we fix all the other parameters of
our model and we train it only once. Thus, we can compare the results for the
different values of the regularization parameter. In Fig. 2, we present the
evolution of the loss function during the training process for the training
and the test set of the Lorenz63 system.

Figure 2: The evolution of the training (red dashed line) and the test set
(solid black line) loss in logarithmic scale, as a function the number of
epochs of training for the Lorenz63 system.

###  III.3 Determining $\alpha$-dependent latent space dimension

The regularization parameter, $\alpha$, controls the size of the latent space
and plays a crucial role in determining the number of the important nodes in
the latent space. In our approach, we define as important the nodes in the
latent space that output a value larger than the regularization parameter,
$\alpha$. Thus, we calculate the average of non-zero nodes in the test set as
a function of $\alpha$. The results for the three systems in consideration are
shown in Fig. 3.

###  III.4 Largest Lyapunov Exponent as a Metric

The Lyapunov exponents are an essential concept of nonlinear dynamics and
present specific features of chaotic systems. Two trajectories of the
attractor initiated at points displaced infinitesimally from each other will
diverge exponentially if the attractor is chaotic. The exponential discrepancy
is described by the largest Lyapunov exponent (LLE). This definition can be
extended to the spectrum of Lyapunov exponents. The LLE of a system can be
calculated using two approaches. The first approach uses a Jacobian matrix and
can be applied when the mathematical model describing the evolution of the
system is known. This method gives the Lyapunov spectrum. The second one can
be used to estimate the LLE of a time series by applying the definition of the
local divergence on a segment of a given trajectory. Estimating the LLE from a
time series faces several numerical difficulties, like dealing with the noise
in the data Wolf ; Rosenstein ; Eckmann . The basic definition of the LLE,
$\lambda_{1}$, can be written as:

$d(t)\approx Ce^{\lambda_{1}t}$ (6)

where, $d(t)$, is an average distance at time, $t$ and $C$ is a constant
(usually normalized by the initial separation of neighbors). Here we use the
methodology proposed by Rosenstein, Collins and De Luca Rosenstein , where we
divide the time, $t$, into steps of $\Delta t$, and assume that at time,
$t=i\times\Delta t$, the distance of the $j^{th}$ pair of nearest neighbors,
$d_{j}(i)$, diverges approximately at a rate given by the LLE:

$d_{j}(i)\approx C_{j}e^{\lambda_{1}i\Delta t},$ (7)

where, $C_{j}$ is the initial distance of the selected pair of points. By
taking the logarithm of Eq. (7) we have that:

$\ln d_{j}(i)\approx\ln C_{j}+\lambda_{1}(i\Delta t).$ (8)

Equation (8), represents a set of $M$, (all the values of index $j$)
approximately parallel lines that can be fitted using least squares to
calculate an averaged value and its standard deviation of the LLE.

In this work, we compare the average value and its standard deviation of the
LLE calculated using the original and the reconstructed trajectories. The
reconstructed dynamics obtained by the autoencoder depend on the parameter in
L1 regularization - $\alpha$. We report these findings and the dependence of
LLE on the embedding dimension window in the results section.

##  IV Results

The basic results of this work consist on the evaluation of the effective size
of the latent space of the chaotic autoencoder constructed for each of the
three systems we considered. The size of the latent space depends on the
regularization parameter $\alpha$ and the window size of input sequence of the
time series, $W$. Once the autoencoder is trained we also compare the largest
Lyapunov value it produces.

###  IV.1 Latent space size as a function of the regularization parameter,
$\alpha$

We apply the methods outlined above on the three chaotic systems, presented in
Table 1. In Fig. 3, we present the average number of non-zero nodes at the
latent space (nodes with output value larger than $\alpha$) as well as the
loss of the reconstructed sequence, with respect to the input sequence, as a
function of the L1 regularization parameter, $\alpha$.

Figure 3: The average number of non-zero nodes at the latent space (navy blue
solid line), including the standard deviation of it (error bars) and the
corresponding value of the loss function (red dashed line) of the test data as
a function of the L1 regularization parameter, $\alpha$, for the three
systems.

We conclude that we find the same patterns for all studied systems, i.e., the
loss function increases with the increment of the regularization parameter.
Depending on which level of accuracy could be needed the number of nodes might
change. Especially, we observe that for the Rössler system the autoencoder
needs two nodes at the latent space to reconstruct the input sequence with
error less than $10^{-4}$, while three to four nodes are needed in the case of
the Lorenz systems.

To further demonstrate the efficiency of the reconstructed trajectories by the
autoencoder model we calculate the largest Lyapunov exponent for the input and
the reconstructed sequences of each one of the regularization parameter
values. We calculate the LLE by creating two trajectories of each attractor.
We choose a starting point for one of them and then we displace it by
$10^{-7}$ for the second trajectory. We create the two trajectories by solving
the equations of each attractor, and then we use the trained autoencoder of
the corresponded system to reconstruct the trajectories. In Table 2, we
present the LLE and its standard deviation in parenthesis for the input
sequence, and the reconstructed sequence of each one of the values of the
regularization parameter, $\alpha$. The results indicate that the
reconstructed sequences created by models with $\alpha$ smaller or equal to
$10^{-4}$ have the same characteristics as the sequences obtained by solving
the deterministic system of equations. The error in the average value of LLE
is less than $\pm$ 2%.

The obtained agreement between the system and chaotic autoencoder LLE is quite
promising since it shows that the latter acts as an effectively linear system
driven by a chaotic signal. It thus, demonstrates further that the chaotic
autoencoder is a faithful reconstruction of the original system in the
specific parameter regime.

|  | Rössler | Lorenz63 | Lorenz96
---|---|---|---|---
Sequence | $\alpha$ | X | Y | Z | X | Y | Z | 19 | 20 | 21
Input | - | 1.73 (0.08) | 1.72 (0.05) | 1.74 (0.09) | 1.17 (0.30) | 1.17 (0.35) | 1.14 (0.18) | 3.02 (0.88) | 3.03 (0.35) | 3.06 (0.57)
Reconstructed | $10^{-3}$ | 1.76 (0.05) | 1.73 (0.02) | 1.77 (0.13) | 1.16 (0.42) | 1.13 (0.21) | 1.14 (0.05) | 2.99 (0.53) | 2.63 (0.73) | 2.49 (0.89)
$10^{-4}$ | 1.74 (0.02) | 1.73 (0.02) | 1.73 (0.01) | 1.14 (0.13) | 1.12 (0.11) | 1.14 (0.34) | 3.08 (0.57) | 2.79 (0.53) | 3.03 (0.57)
$10^{-5}$ | 1.77 (0.09) | 1.73 (0.06) | 1.74 (0.09) | 1.14 (0.16) | 1.17 (0.34) | 1.13 (0.15) | 2.88 (0.57) | 2.83 (0.44) | 3.01 (0.54)
$10^{-6}$ | 1.78 (0.09) | 1.73 (0.04) | 1.78 (0.01) | 1.14 (0.26) | 1.17 (0.36) | 1.14 (0.24) | 2.88 (0.63) | 3.01 (0.22) | 2.90 (0.42)
$10^{-7}$ | 1.76 (0.07) | 1.74 (0.02) | 1.78 (0.05) | 1.15 (0.29) | 1.17 (0.26) | 1.13 (0.35) | 2.97 (0.41) | 2.83 (0.56) | 3.04 (0.43)
$10^{-8}$ | 1.73 (0.25) | 1.73 (0.11) | 1.79 (0.05) | 1.17 (0.27) | 1.17 (0.29) | 1.14 (0.07) | 2.87 (0.27) | 3.04 (0.64) | 3.01 (0.50)
Table 2: The largest Lyapunov exponent and its standard deviation in
parenthesis, for the input sequence and for each one of the reconstructed ones
as a function of the regularization parameter, $\alpha$, for the three
systems.

###  IV.2 Latent space dimension as a function of the size of the input
sequence, $W$

It is important to know how the size of the information passed to the
autoencoder (size of input sequence) affects the performance of the model. In
Fig. 4, we present the average number of non-zero nodes in the latent space
that were used to reconstruct the input sequence as a function of the size of
the input sequence, $W$. We fix the regularization parameter $\alpha$ to
$10^{-5}$ and range the size of the input sequence, $W$ from 9 to 37 time
steps using steps of 7 time steps.

Figure 4: The average number of non-zero nodes at the latent space (navy blue
solid line), including the standard deviation of it (error bars) and the
corresponding value of the loss function (red dashed line) of the test data as
a function of the size of the input sequence, $W$, for the three systems.

In the case of the Rössler system, the autoencoder needs two non-zeros nodes
in the latent space to reconstruct a sequence of more than 16 time steps.
Three non-zero nodes are needed in the Lorenz63 model to reconstruct a
sequence of 23 time steps. We observe again, that the Lorenz96 system is more
complex than the other two, and the adequate number of nodes to represent it
depends strongly on the desired accuracy. In Table 3, we present the LLE and
its standard deviation in parenthesis for the input sequence, and the
reconstructed sequence of each one of the values of the size of the input
sequence, $W$. The error in the average value of LLE found to be less than 5%
in all cases.

|  | Rössler | Lorenz63 | Lorenz96
---|---|---|---|---
Sequence | $W$ | X | Y | Z | X | Y | Z | 19 | 20 | 21
Input | - | 1.73 (0.08) | 1.72 (0.05) | 1.74 (0.09) | 1.17 (0.30) | 1.17 (0.35) | 1.14 (0.18) | 3.02 (0.88) | 3.03 (0.35) | 3.06 (0.57)
Reconstructed | 9 | 1.77 (0.02) | 1.75 (0.04) | 1.70 (0.19) | 1.15 (0.37) | 1.17 (0.34) | 1.15 (0.16) | 2.92 (0.40) | 2.92 (0.64) | 3.95 (0.51)
16 | 1.73 (0.21) | 1.72 (0.34) | 1.80 (0.12) | 1.13 (0.18) | 1.13 (0.34) | 1.13 (0.40) | 2.95 (0.64) | 3.05 (0.57) | 3.08 (0.21)
23 | 1.73 (0.21) | 1.76 (0.08) | 1.67 (0.44) | 1.14 (0.40) | 1.15 (0.32) | 1.15 (0.51) | 2.98 (0.39) | 2.92 (0.59) | 3.08 (0.38)
30 | 1.76 (0.14) | 1.73 (0.05) | 1.73 (0.85) | 1.14 (0.29) | 1.18 (0.13) | 1.18 (0.44) | 2.94 (0.52) | 2.97 (0.76) | 2.92 (0.41)
37 | 1.74 (0.16) | 1.74 (0.13) | 1.82 (0.79) | 1.14 (0.18) | 1.16 (0.44) | 1.17 (0.41) | 2.96 (0.35) | 3.03 (0.36) | 3.02 (0.46)
Table 3: The largest Lyapunov exponent and its standard deviation in
parenthesis, for the input sequence and for each one of the reconstructed ones
as a function of the size of the input sequence, $W$, for the three systems.

##  V Conclusions

In this work, we attempted to connect chaotic dynamical systems with neural
network-based autoencoders and find the minimal set of NN variables that may
reproduce these dynamical systems. The construction of chaotic autoencoders is
performed through deep learning networks, and in the process of their use, an
information bottleneck is formed. The latter is determined exclusively by the
dimension of the latent space. We aimed to find the minimal dimension of the
latent space that faithfully reproduces the chaotic systems we studied. This
dimension depends on the required accuracy at the output. Furthermore, it also
depends on the window size we choose in order to perform the attractor
reconstruction. The analysis of the latent space dimensionality as a function
of the regularization parameter and the input sequence size gives compatible
numbers and this feature adds to the consistency of the picture. In addition,
the resulting maximal Lyapunov exponent of the chaotic autoencoder is
respectively very close to each of the chaotic systems we analyzed. This is a
significant outcome since it shows that the chaotic autoencoder can represent
faithfully complex, irregular dynamics.

An important issue that opens up from the present work is the connection of
the dimensionallity of the latent space to that of the embedding dimension of
the time series. Intuitively one might think that the latent space is in some
sense the minimal space required to represent a chaotic dynamical system. As a
result, the latent space dimensions should be smaller or at least equal to the
embedding dimension since the latter gives the framework space for analyzing
the chaotic series. We observe that this feature is true in the examples we
addressed since we find latent space dimension smaller or equal to the number
of differential equations that generate respective data. An additional
appealing outcome of the present analysis is that the largest Lyapunov
exponent that is connected with the autoencoder system appears to be very
close to the one determined directly from the dynamical system. Clearly, these
two exponents are evaluated differently since the chaotic autoencoder cannot
by construction generate its dynamics and depends on the input time series.
The Lyapunov exponent, in this case, is obtained by observing the proximity of
two nearby chaotic series with a very small difference in their initial
conditions. However, the autoencoder can generate similar divergence in the
difference of the two series as that of the original system shows that it
provides a faithful representation of the original system. In other words, we
find that the chaotic autoencoder provides a faithful representation of time
series segments and the more intricate internal dynamics of the chaotic
system.

The present work opens up an information window in dynamic systems. It
demonstrates that AI methods, when applied properly, may give information on
complex dynamical systems and their various reduced representations. This
information is both quantitative and provides additional representations of
dynamical systems. In some sense, it also shows that neural networks are
ubiquitous and may capture accurately a high degree of information contained
in complex dynamical systems. This result treated a bit more philosophically
may be seen as a rather optimistic one for us since it demonstrates that human
neural networks are indeed capable of capturing the complexity of the outside
world.

###### Acknowledgements.

We acknowledge the co-financing of this research by the European Union and
Greek national funds through the Operational Program Crete 2020-2024, under
the call "Partnerships of Companies with Institutions for Research and
Transfer of Knowledge in the Thematic Priorities of RIS3Crete", with project
title “Analyzing urban dynamics through monitoring the city magnetic
environment” (project KPHP1 - 0029067) and also by the Ministry of Science and
Higher Education of the Russian Federation in the framework of Increase
Competitiveness Program of NUST "MISiS" (No. K2-2019-010), implemented by a
governmental decree dated 16th of March 2013, N 211.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * (1) V. S. Anishchenko, V. Astakhov, A. Neiman, T. Vadivasova, L. S.-Geier, Nonlinear Dynamics of Chaotic and Stochastic Systems, Springer-Verlag Berlin Heidelberg, $2^{nd}$ edition, 0172-7389, (2007).
  * (2) H. Abarbanel, Analysis of Observed Chaotic Data, Springer-Verlag New York, 1996.
  * (3) H. Kantz and T. Schreiber, Nonlinear Time Series Analysis, 7, Cambridge University Press, (2004).
  * (4) E. Bradley and H. Kantz, "Nonlinear time-series analysis revisited," Chaos, 25, 097610, (2015).
  * (5) H. Abarbanel, R. Brown, J. J. Sidorowich and L. Sh. Tsimring, Rev. Mod. Phys. 65, 1331, (1993).
  * (6) E. R. Deyle, G. Sugihara, Generalized Theorems for Nonlinear State Space Reconstruction, PLOS ONE, 6, (3), (2011)
  * (7) P. Verdes, P. Granitto, H. Navone and H. Ceccatto, "Forecasting chaotic time series: Global vs. local methods," Novel Intelligent Automation and Control Systems, 1, 129, (1998).
  * (8) I. Goodfellow, Y. Bengio, A. Courville and Y. Bengio, Deep Learning, MIT Press, 1, (2016).
  * (9) M. I. Jordan and T. M. Mitchell, "Machine learning: Trends, perspectives, and prospects", Science, 349, 6245, pp. 255, (2015).
  * (10) A. Wolf, J.B. Swift, H.L. Swinney and J.A. Vastano, "Determining Lyapunov exponents from a time series," Physica D, 16, p. 285, (1985).
  * (11) M. Sano and Y. Sawada, "Measurement of the Lyapunov spectrum from a chaotic time series," Phys. Rev. Lett. , 55, p. 1082 , (1985).
  * (12) M. T. Rosenstein, J. J. Collins and C. J. De Luca, "A practical method for calculating largest Lyapunov exponents," Physica D: Nonlinear Phenomena, 65, (1-2), p. 117–134, (1993).
  * (13) J. P. Eckmann, S. O. Kamphorst, D. Ruelle, and S. Ciliberto, "Liapunov Exponents from Time Series," Phys. Rev. A 34, pp. 4971-4979 , (1986).
  * (14) U. Parlitz, "Identification of true and spurious Lyapunov exponents from time series," International journal of bifurcation and chaos, 01, 155, (1992).
  * (15) M. Nielsen, Neural Networks and Deep Learning, San Francisco: Determination Press, (2015).
  * (16) Lu Z, Pathak J, Hunt B, Girvan M, Brockett R and Ott E. , Reservoir observers: model-free inference of unmeasured variables in chaotic systems, Chaos, 27, 041102, (2017).
  * (17) Budhiraja R., Kumar M., Das M. K., Bafila A. S. , Singh S. , "A reservoir computing approach for forecasting and regenerating both dynamical and time-delay controlled financial system behavior.," PLoS ONE, 16, no. 2, (2021).
  * (18) P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis and P. Koumoutsakos, "Data-Driven Forecasting of High-Dimensional Chaotic Systems with Long Short-Term Memory Networks," Proc R Soc A. 474 , p. 2213, (2018).
  * (19) G. Neofotistos, M. Mattheakis, G.D. Barmparis, J. Hizanidis, G.P. Tsironis, E. Kaxiras, Machine learning with observers predicts complex spatiotemporal behavior, Front. Phys. 7, 24, (2019).
  * (20) S. L. Brunton, J. L. Proctor, J. N. Kutz, "Discovering governing equations from data by sparse identification of nonlinear dynamical systems," Proc. Natl. Acad. Sci. U.S.A., 113, p. 3932–3937, (2016).
  * (21) F.Dercole, M.Sangiorgio and Y.Schmirander, "An empirical assessment of the universality of ANNs to predict oscillatory time series," IFAC-PapersOnLine, 53, 2, pp. 1255-1260, (2020).
  * (22) B. Lusch, J. N. Kutz, S. L. Brunton, " Deep learning for universal linear embeddings of nonlinear dynamics," Nat. Commun, 9, p. 4950, (2018).
  * (23) Lopez Pinaya W. H., Vieira S., Garcia-Dia, R. and Mechelli A., "Autoencoders. Machine Learning," p. 193–208, (2020).
  * (24) M. Kramer, "Nonlinear principal component analysis using autoassociative neural networks," AIChE Journal, 37, 2, p. 233–243, (1991).
  * (25) T. Sauer, J. A. Yorke and M. Casdagli, "Embedology," Journal of Statistical Physics, 65 , 3, p. 579–616, (1991).
  * (26) R. Hilborn, Chaos and nonlinear dynamics: an introduction for scientists and engineers, Oxford University Press, (2000).
  * (27) E. N. Lorenz, "The Local Structure of a Chaotic Attractor in Four Dimensions," Physica D, 13, pp. 90-104 , (1984).
  * (28) D. S. Broomhead and G. P. King, "Extracting qualitative dynamics from experimental data," Physica D: Nonlinear Phenomena, 20, 2, pp. 217-236, (1986).
  * (29) F. Takens, "Dynamical Systems and Turbulence," Lecture Notes in Mathematics (D. A. Rand and L. S. Young, eds), Berlin, Springer Verlag, pp. 366-381, (1981).
  * (30) N. H. Packard, J. P. Crutchfield, J. D. Farmer and R. S. Shaw, "Geometry from a Time Series," Phys. Rev. Lett. 45, (9), p. 712–716, (1980).
  * (31) M. T. Rosenstein, J. J. Collins and C. J. De Luca, Reconstruction expansion as a geometry-based framework for choosing proper delay times, Phys. D 73, 1–2, 82–98, (1994).
  * (32) D.P. Kingma and J. Ba, Adam: a method for stochastic optimization, arXiv:1412. 6980v9, preprint, (2017).
  * (33) Xavier Glorot and Yoshua Bengio, Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, PMLR 9:249-256, (2010).
  * (34) Abadi, Martín and Barham, Paul and Chen, Jianmin and Chen, Zhifeng and Davis, Andy and Dean, Jeffrey and Devin, Matthieu and Ghemawat, Sanjay and Irving, Geoffrey and Isard, Michael and others, Tensorflow: A system for large-scale machine learning, 12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16), 265, (2016).
  * (35) Chollet, François et al., Keras, https://keras.io, (2015)


# Automaton of molecular perceptions in biochemical reactions

Stefano Maestri<EMAIL_ADDRESS>Emanuela Merelli  School of Science
and Technology, University of Camerino, Camerino, 62032, Italy

###### Abstract

Local interactions among biomolecules, and the role played by their
environment, have gained increasing attention in modelling biochemical
reactions. By defining the automaton of molecular perceptions, we explore an
agent-based representation of the behaviour of biomolecules in living cells.
Our approach considers the capability of a molecule to perceive its
surroundings a key property of bimolecular interactions, which we investigate
from a theoretical perspective. Graph-based reaction systems are then
leveraged to abstract enzyme regulation as a result of the influence exerted
by the environment on a catalysed reaction. By combining these methods, we aim
at overcoming some limitations of current kinetic models, which do not take
into account local molecular interactions and the way they are affected by the
reaction environment.

###### keywords:

agent-based modelling, graph-based reaction systems , bimolecular interactions
, finite-state automata

##  1 Introduction

Biochemical reactions, especially those involved in metabolic and signaling
pathways, have been widely studied through computational approaches [1, 2, 3,
4]. In particular, Systems Biology provides computational models focused
mostly on the kinetic properties of the enzymatic reactions, represented
through sets of differential equations [5]; this approach may constitute a
limitation, especially since kinetic parameters sampled in vitro might not
capture the dynamic interactions that enzymes carry out with their environment
[6]. In this context, the latter is also considered as homogeneous, and
enzymes as functions from reactants to products, overlooking the granularity
that characterises molecular interactions.

In recent years, this problem has been tackled by modelling and simulating the
molecules involved in a biochemical reaction as agents, autonomous systems
able to perceive and interact with the other agents populating their
environment [7, 8]. Such an approach allows capturing the fundamental role
that molecular perception has on the efficiency of a biochemical reaction, a
property that has been experimentally observed in the form synchronised
oscillations among interacting molecules [9].

According to the Michaelis-Menten model of enzyme kinetics, an enzymatic
reaction can be represented as:

$E+S\leftrightharpoons ES\rightarrow E+P$ (1)

where E is an enzyme, S its substrate and P the product of the reaction
catalysed by E; assuming the steady-state approximation, we can consider ES as
constant [10, 11].

As part of an agent-based model of a biochemical reaction, we described this
pattern through a basic reaction automaton [12]. It was intended just to
support the interaction phases accounted in the agent-based model and does not
explicitly consider agent’s perception, as well as ignores the effects of
enzyme regulation on the generation of the reaction products.

In this manuscript, we refine this idea from a theoretical perspective, by
considering the environment as a first-class component of the system, which is
perceived by the agents representing the molecule involved in the reaction; we
also leverage the capabilities of the graph-based reaction systems [13] to
move a first step in modelling the regulation carried out, on the enzymatic
activity, by specific molecules of the environment.

##  2 Materials and Methods

The work described in this manuscript relies on automata theory, agent-based
modelling and graph-based reaction systems. In this section, we provide a
brief introduction to graph-based reaction systems, which, among the other, is
a relatively new theory, and maybe less known.

Reaction systems were introduced as a formal framework to model biochemical
reactions in living cells. Although the original approach is purely set-
theoretic [14], a graph-based description is provided by Kreowski and
Rozenberg [13], allowing to deal with graph-related problems, such as those we
tackle in this manuscript. We report here some basic concepts, as described in
the aforementioned source article, needed for the theoretical constructions
proposed in the next section.

According to Kreowski and Rozenberg [13]:

  * •

A directed and edge-labelled graph is a system $G=(V,\Sigma,E)$, where $V$ is
a finite set of nodes, $\Sigma$ is a finite set of edge labels, and
$E\subseteq V\times V\times\Sigma$ is a set of edges.

  * •

Given the graphs $H$ and $G$, such that $\Sigma_{H}=\Sigma_{G}$, $H$ is a
subgraph of $G$ if $V_{H}\subseteq V_{G}$, and $E_{H}\subseteq E_{G}$. We
denote with $Sub(G)$ the set of all the subgraphs of $G$.

  * •

A selector $\mathcal{S}$ of a graph $G=(V,E)$, is an ordered pair
$\mathcal{S}=(X,Y)$ such that $X\subseteq V$ and $Y\subseteq E$.

  * •

The extraction $U(H)$ of $H\in Sub(G)$ is the pair $U(H)=(V_{H},E_{H})$,
respectively of the set of nodes and set of edges of the graph $H$.

  * •

A reaction b over a non-empty, directed, and edge-labelled graph
$B=(V_{B},\Sigma_{B},E_{B})$, called the background graph, is a triple
$b=(R,I,P)$, where $R$ and $P$ are non-empty subgraphs of $B$, and $I$ is a
selector of $B$ such that $I\cap\,U(R)=(\emptyset,\emptyset)$.

  * •

Let $T$ be a subgraph of $B$.

    * –

A reaction $b=(R,I,P)$ over $B$ is enabled by $T$, denoted by $en_{b}(T)$, if
$R\in Sub(T)$ and $I\cap U(T)=(\emptyset,\emptyset)$.

    * –

The result of a reaction $b$ on $T$, denoted by $res_{b}(T)$, is defined by
$res_{b}(T)=P_{b}$ if $en_{b}(T)$; $res_{b}(T)=\emptyset$ otherwise.

    * –

The result of a set of reactions $A$ over $B$ on $T$ , denoted by
$res_{A}(T)$, is defined by: $res_{A}(T)=\bigcup_{b\in A}res_{b}(T)$.

    * –

A graph-based reaction system is a pair $\mathcal{A}=(B,A)$, where $B$ is a
finite non-empty graph, and $A$ is a set of reactions over $B$. Also,
$\mathcal{A}$ induces the function $res_{\mathcal{A}}:Sub(B)\rightarrow
Sub(B)$, called the result function of $\mathcal{A}$, such that, for each
state $T\in Sub(B)$, $res_{\mathcal{A}}(T)=res_{A}(T)$.

##  3 Results

###  3.1 A basic enzymatic reaction model

An enzymatic reaction, as introduced in Section 1, can be modelled by
explicitly taking into account the possible role of coenzymes, such as ATP or
NADH.

Let $\Sigma=\\{c,p,r,s\\}$ be an alphabet, where $c$ represents a coenzyme,
$p$ the main product of the reaction, $r$ the result of the conversion of a
coenzyme (if any), and $s$ a substrate molecule. We clarify that a coenzyme
is, indeed, substrate of the related enzyme, but we use the above naming to
emphasise its specific role of supporting the reaction. Let also
$Q=\\{E,Ec,Ep,Er,Es,Esc\\}$ be a set of states, where $E$ represents a free
enzyme, $Ec$ the binding of a coenzyme, $Es$ the enzyme bound to a substrate
molecule, $Esc$ an enzyme saturated by a substrate molecule and a coenzyme;
$Ep$ and $Er$ are the states through which, respectively, the main product of
the reaction and the product of the coenzyme conversion are released. Binding
a substrate molecule and a coenzyme can happen without a specific order,
however there is no case in which a coenzyme can directly generate a product.
We also assume that, when the reaction has two products, they are not released
exactly in the same instant, even though, from the biological perspective,
this time delay may be negligible. The basic behaviour of an enzymatic
reaction can thus be described by the deterministic finite state automaton
$\mathcal{R}=(Q,\Sigma,\delta,s_{0},\\{f\\})$, with $s_{0}=f=E$ being the
initial and final state, and $\delta:Q\times\Sigma\rightarrow Q$ the state
transition relation defined as follows:

State Input | $\bm{s}$ | $\bm{c}$ | $\bm{p}$ | $\bm{r}$
---|---|---|---|---
$\bm{E}$ | $Es$ | $Ec$ | $-$ | $-$
$\bm{Es}$ | $-$ | $Esc$ | $E$ | $-$
$\bm{Ec}$ | $Esc$ | $-$ | $-$ | $-$
$\bm{Esc}$ | $-$ | $-$ | $Er$ | $Ep$
$\bm{Ep}$ | $-$ | $-$ | $E$ | $-$
$\bm{Er}$ | $-$ | $-$ | $-$ | $E$

Note that, from $Es$ (i.e., the state in which the enzyme $E$ is bound to the
substrate molecule s), the automaton can reach either the state $E$, by
generating the product $p$, or the state $Esc$, by binding the coenzyme $c$.
This approach allows the enzymatic reaction automaton to generalise the
behaviour of enzymes that require a coenzyme to catalyse the reaction, as well
as of those able to directly facilitate the substrate-to-product conversion.

The automaton state graph is shown in Fig. 1.

$E$$Ec$$Es$$Esc$$Ep$$Er$$c$$s$$s$$c$$r$$p$$p$$r$$p$ Figure 1: Automaton
representing the basic behaviour of an enzymatic reaction, which may or may
not involve a coenzyme.

This model of enzymatic reaction is a refined version of the one described in
[15], and does not consider the effect of molecular perception. However, form
an agent-based perspective, enzymes are active entities able to perceive and
modify their environment, which is populated by agents modelling substrate
molecules and other enzymes.

###  3.2 Enhancing the enzymatic reaction model with perception

To endow $\mathcal{R}$ with perception, we should add a level of
nondeterminism, and introduce $\epsilon$-transitions. The latter are needed to
model the unstable states in which an enzyme perceives a cognate molecule, but
the actual binding is not yet certain.

We therefore construct the extended reaction automaton
$\mathcal{R^{{}^{\prime}}}$, which differs from the previous one in its set of
states and transition function.

######  Definition 3.1 (Perception-based reaction automaton).

We define perception-based reaction automaton the 5-tuple
$\mathcal{R^{{}^{\prime}}}=(Q^{{}^{\prime}},\Sigma\cup\\{\epsilon\\},\delta^{{}^{\prime}},s_{0},\\{f\\})$,
such that:

  * •

$Q^{{}^{\prime}}=\mathscr{S}U\mathscr{P}$, where
$\mathscr{S}=\\{E,Ec,Ep,Er,Es,ES\\}$ and
$\mathscr{P}=\\{E_{c},Ec_{s},E_{s},Es_{c}\\}$.

    * –

$\mathscr{S}$ is the set of stable states, where an enzyme $E$ is free, bound
to a cognate molecule or coenzyme in a molecular complex ($Es$, $Ec$), or is
saturated ($ES$) and therefore ready to catalyse the reaction, release its
product (through the $Er$ and $Ep$ states, when a coenzyme is involved), and
thus return free;

    * –

$\mathscr{P}$ represents the set of perceiving states, where the free enzyme
(or a complex) perceives a cognate molecule $s$ or coenzyme $c$ (indicated as
subscripts).

  * •

$\Sigma=\\{c,p,r,s\\}$ and $s_{0}=f=E$

  * •

$\delta^{{}^{\prime}}:Q\times\Sigma\cup\\{\epsilon\\}\rightarrow 2^{Q}$, is
given by the following transition table:

State Input | $\bm{s}$ | $\bm{c}$ | $\bm{p}$ | $\bm{r}$ | $\bm{\epsilon}$
---|---|---|---|---|---
$\bm{E}$ | $\\{E_{s}\\}$ | $\\{E_{c}\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{E\\}$
$\bm{Es}$ | $\\{\\}$ | $\\{Es_{c}\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{Es\\}$
$\bm{Ec}$ | $\\{Ec_{s}\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{Ec\\}$
$\bm{ES}$ | $\\{\\}$ | $\\{\\}$ | $\\{Er,E\\}$ | $\\{Ep\\}$ | $\\{\\}$
$\bm{Ep}$ | $\\{\\}$ | $\\{\\}$ | $\\{E\\}$ | $\\{\\}$ | $\\{\\}$
$\bm{Er}$ | $\\{\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{E\\}$ | $\\{\\}$
$\bm{E_{s}}$ | $\\{Es,ES\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{E\\}$
$\bm{Es_{c}}$ | $\\{\\}$ | $\\{ES\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{Es\\}$
$\bm{E_{c}}$ | $\\{\\}$ | $\\{Ec\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{E\\}$
$\bm{Ec_{s}}$ | $\\{ES\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{\\}$ | $\\{Ec\\}$

As mentioned above, the $\epsilon$-transition allows a stable enzyme or
complex to not perceive any cognate molecule in its surrounding environment or
to return, if the binding of the perceived molecule does not happen, from a
perceiving state to a previous stable state.

Introducing perceiving states in the model forces the nondeterminism on the
pairs $(E_{s},s)$ and $(ES,p)$ to distinguish the case in which the reaction
involves a coenzyme from that where the sole substrate molecule $s$ is
converted.

The $ES$ state represents a saturated enzyme, which is not able to perceive
its surroundings; similarly, the enzyme has no perception in the $Er$ and $Ep$
states, where it releases the products of the reaction. This behaviour is
represented by the state graph in Fig. 2.

$E$$E_{s}$$E_{c}$$Ec$$Es$$Ec_{s}$$Es_{c}$$ES$$Er$$Ep$$s$$c$$\epsilon$$\epsilon$$c$$s$$\epsilon$$\epsilon$$s$$c$$\epsilon$$\epsilon$$c$$s$$p$$r$$r$$p$$s$$p$$\epsilon$

Figure 2: State graph of the perception-based reaction automaton; it
represents the generalised behaviour of the perception-based enzymatic
reactions (see Definition 3.3).

######  Definition 3.2 (Enzyme perception and reaction anabler).

Given a perception-based reaction automaton
$\mathcal{R^{{}^{\prime}}}=(Q^{{}^{\prime}},\Sigma\cup\\{\epsilon\\},\delta^{{}^{\prime}},s_{0},\\{f\\})$,
and a set $\Pi\subset\Sigma/\\{p\\}$, the function
$\mu:\mathscr{S}\times\Pi\rightarrow\mathscr{P}$ is called enzyme perception
in $\mathcal{R^{{}^{\prime}}}$. The set $\Pi$, defined reaction enabler,
contains all the cognate molecules that the enzyme must perceive in its
surrounding environment in order to proceed during the reaction.

It is possible to identify an equivalence relation $\equiv$ that characterises
the equivalence class $[\mathcal{R^{{}^{\prime}}}]_{\equiv}$ of the automaton
$\mathcal{R^{{}^{\prime}}}$. Although a rigorous definition of this class is
beyond the scope of the present article (and will be provided in a future
work), intuitively, we can say that it represents the class of all the
automata that show the same behaviour of $\mathcal{R^{{}^{\prime}}}$.

######  Definition 3.3 (Perception-based enzymatic reaction).

Being $[\mathcal{R^{{}^{\prime}}}]_{\equiv}$ the equivalence class of the
perception-based automaton $\mathcal{R^{{}^{\prime}}}$, a 5-tuple
$\mathcal{R}_{p}=(Q,\Sigma,\delta,s_{0},F)$ is a perception-based enzymatic
reaction iff $\mathcal{R}_{p}\in[\mathcal{R^{{}^{\prime}}}]_{\equiv}$ (that
is, $\mathcal{R}_{p}\equiv\mathcal{R^{{}^{\prime}}}$).

######  Example 3.1.

The phosphorylation of fructose 6-hosphate (F6P) to fructose 1,6-bisphosphate
(F16bP), which is carried out by the phosphofructokinase (PFK) and coupled
with the hydrolysis of ATP to ADP, is modelled through the perception-based
enzymatic reaction shown in Fig. 3.

$PFK$$PFK_{F6P}$$PFK_{ATP}$$PFK$$ATP$$PFK$$F6P$$PFK$$ATP$F6P$PFK$$F6P$ATP$PFK$$ATPF6P$$PFK$$ADP$$PFK$$F16bP$$F6P$$ATP$$\epsilon$$\epsilon$$ATP$$F6P$$\epsilon$$\epsilon$$F6P$$ATP$$\epsilon$$\epsilon$$ATP$$F6P$$F16bP$$ADP$$ADP$$F16bP$$\epsilon$

Figure 3: Perception-based enzymatic reaction performed by phosphofructokinase
(PFK).

###  3.3 Enzyme regulation in perception-based reactions

An important property of the enzymatic reactions, especially those involved in
metabolic and signaling pathways, is the capability of being regulated (that
is, activated or inhibited) in the process of generating their products.
Enzyme regulation is often significantly complex; we propose here a first
approach to tackle this property in modelling a perception-based enzymatic
reaction. The core idea is to abstract the inhibition process as a “switch”
that completely disable the reaction, so that it cannot generate its products;
to implement this idea, we rely on the framework of graph-based reaction
systems, treated in [13], and outlined in Section 2 of this manuscript.

According to [13], a reaction over a non-empty, directed and edge-labelled,
graph $B=(V_{B},\Sigma_{B},E_{B})$, called background graph, is a triple
$b=(R,I,P)$, where $R$ and $P$ are non-empty subgraphs of $B$, and
$I=(V_{I},E_{I})$ is a selector of $B$ such that
$I\cap\,U(R)=(\emptyset,\emptyset)$. $R$ is the reactant graph of $b$, $I$ is
the inhibitor of $b$, and $P$ is the product graph of $b$.

We recall that, a selector $\mathcal{S}$ of a graph $G=(V,E)$, is an ordered
pair $\mathcal{S}=(X,Y)$ such that $X\subseteq V$ and $Y\subseteq E$, and that
the extraction $U(H)$ of $H\in Sub(G)$ is the pair $U(H)=(V_{H},E_{H})$,
respectively of the set of nodes and set of edges of the graph $H$.

In our setting, the background graph $B$ represents all the possible spatial
configurations that the molecules can assume in the modelled portion of
cytoplasm. Each of them constitutes the environment of the perception-based
enzymatic reactions.

To introduce enzyme perception in a graph-based reaction, we need to define a
new selector $\mathcal{S}_{\Pi}$ (perception selector), which identifies all
the edges of $B$ labelled with the substrate molecules enabling a perception-
based enzymatic reaction $\mathcal{R}_{b}$.

######  Definition 3.4 (Perception selector).

Being $\mathcal{R}_{p}=(Q,\Sigma,\delta,s_{0},F)$ a perception-based enzymatic
reaction, $B=(V_{B},\Sigma_{B},E_{B})$ a background graph, and
$\Upsilon=\\{I\,|\,I=(V_{I},E_{I})\in b,\forall$ $b=(R,I,P)$ over $B\\}$, such
that $\Sigma_{B}=\Sigma\cup\Gamma$, where $\Gamma=\\{i\,|\,(v,v^{\prime},i)\in
E_{I}$ and $v,v^{\prime}\in V_{I},\forall\,I=(V_{I},E_{I})\in\Upsilon\\}$.
Being also $\Pi\subset\Sigma$ the reaction enabler of $\mathcal{R}_{b}$, and
$E_{\Pi}\subset E_{B}$ such that $E_{\Pi}=\\{(v,v^{\prime},x)\,|\,x\in\Pi\\}$,
a perception selector of $\mathcal{R}_{p}$ is a pair
$\mathcal{S}_{\Pi}=(\\{v,v^{\prime}\,|\,(v,v^{\prime},x)\in
E_{\Pi}\\},E_{\Pi}).$

Now, we can extend the enabled reactions defined in [13] (see Section 2) with
the concept of enzyme perception.

######  Definition 3.5 (Perception-enabled reaction, result by perception).

Given a graph $T\in Sub(B)$, a reaction $b=(R,I,P)$ over $B$ is perception-
enabled by $T$, denoted $pen_{b}(T)$, iff

  * •

$R,P\in Sub(T)$ and $U(R)\cap\mathcal{S}_{\Pi}\neq(\emptyset,\emptyset)$;

  * •

$U(T)\cap I=(\emptyset,\emptyset)$.

Otherwise, b is disabled by $T$, denoted $dis_{b}(T)$.

The result by perception $resp_{b}(T)$ of a reaction $b$ on $T$, is defined as
$resp_{b}(T)=P_{b}$ iff $pen_{b}(T)$; if $dis_{b}(T)$, then
$resp_{b}(T)=\emptyset_{\Sigma}$. The result by perception of a set of
reactions $A$ over $B$ on $T$, denoted by $resp_{A}(T)$, is defined as
$resp_{A}(T)=\bigcup_{b\in A}resp_{b}(T)$.

The subset $T$ of $B$ represents the configuration that the molecules in the
modelled cytoplasm portion assume at a specific time (which may be, as an
example, the time-step of an agent-based simulation). A perception-based
reaction $\mathcal{R}_{b}$ is enabled by this environment configuration only
if the enzyme that carries out $\mathcal{R}_{b}$ perceives its cognate
molecules (i.e., the substrate molecules that enable the reaction) and if it
is not inhibited. In Definition 3.5, such conditions are guaranteed ed by
imposing, respectively, $U(R)\cap\mathcal{S}_{\Pi}\neq(\emptyset,\emptyset)$
and $U(T)\cap I=(\emptyset,\emptyset)$.

######  Remark 3.5.1.

Since $R,P\in Sub(T)\in Sub(B)$ and $\mathcal{S}_{\Pi}$ is defined over
$\Pi\subset\Sigma\subset\Sigma_{B}$, perception-enabled reactions are special
cases of graph-based reactions.

Kreowski and Rozenberg [13] define a graph-based reaction system as a pair
$\mathcal{A}=(B,A)$, where $B$ is a finite non-empty graph, and $A$ is a set
of reactions over $B$.

######  Definition 3.6 (Result by perception function).

A function $resp_{\mathcal{A}}:Sub(B)\rightarrow Sub(B)$, called the result by
perception function of $\mathcal{A}=(B,A)$, is such that, for each state $T\in
Sub(B)$, $resp_{\mathcal{A}}(T)=resp_{A}(T)$.

######  Definition 3.7 (Perception-based enzymatic reaction inhibition).

A perception-based enzymatic reaction
$\mathcal{R}_{p}=(Q,\Sigma,\delta,s_{0},F)$ is inhibited by a configuration of
its environment, defined through a graph T, iff the graph-based reaction
system $\mathcal{A}(\mathcal{R}_{p})=(B(\mathcal{R}_{p}),A(\mathcal{R}_{p}))$,
associated to $\mathcal{R}_{p}$ through the steps described in Section 5 of
[13], is such that $T\in Sub((B(\mathcal{R}_{p}))$, and, $\forall$ reaction
$b=(R,I,P)\in A(\mathcal{R}_{p})$, $\exists I=(V_{I},E_{I})$ such that
$U(T)\cap I\neq(\emptyset,\emptyset)$, that is,
$resp_{b}(T)=\emptyset_{\Sigma}$.

##  4 Conclusions

This manuscript provides a brief introduction on how finite-state automata and
graph-based reaction systems can be combined to model the capabilities of
biomolecules to perceive and interact with their environment.

We start by refining a previously introduced reaction automaton [15, 12], and
by endowing this model with environment perception. Providing molecules with
the ability to perceive cognate partners in their surroundings allows us to
identify the equivalence class of the perception-based reaction automata, to
which we associate graph-based reaction systems [13]. From this new
perspective, the environment is seen as a background graph $B$, while a
molecules’ space configuration (e.g., at a given time) as a subgraph $T$ that
determines if a reaction can be carried out or not according to that specific
context. Given these premises, we are then able to formally define the
enzymatic reaction inhibition, carried out by the reaction environment, as a
consequence of molecular perception.

Our analysis aims to set the basis for a novel modelling approach that, laying
on solid formal frameworks, takes into account the fundamental effects of the
reaction environment on biomolecular interactions.

Although the proposed approach is strongly theoretical, it can be expanded in
order to specify the behaviour of interacting molecules in a multiagent
simulation [12]. Further improvements of this work will also model the
reversibility of perception-enabled reaction systems and introduce a way to
influence the reaction direction through measurement values [16, 17].

## Acknowledgments

This article is the result of the research project funded by the Future and
Emerging Technologies (FET) programme within the Seventh Framework Programme
(FP7) for Research of the European Commission, under the FET-Proactive grant
agreement TOPDRIM (www.topdrim.eu), number FP7-ICT- 318121.

## Author contributions statement

EM designed and supervised the research. SM implemented the method, performed
the research and wrote the paper. All authors revised the paper and approved
the final manuscript.

## Competing interests

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## References

  * [1] H. Berry, Monte Carlo Simulations of Enzyme Reactions in Two Dimensions: Fractal Kinetics and Spatial Segregation, Biophysical Journal 83 (4) (2002) 1891–1901 (Oct. 2002).  doi:10.1016/S0006-3495(02)73953-2.
  * [2] K. Takahashi, S. N. V. Arjunan, M. Tomita, Space in systems biology of signaling pathways - towards intracellular molecular crowding in silico, FEBS Lett. 579 (8) (2005) 1783–1788 (Mar. 2005).  doi:10.1016/j.febslet.2005.01.072.
  * [3] E. Bartocci, F. Corradini, E. Merelli, L. Tesei, Detecting synchronisation of biological oscillators by model checking, Theoretical Computer Science 411 (20) (2010) 1999–2018 (Apr. 2010).  doi:10.1016/j.tcs.2009.12.019.
  * [4] F. Buti, D. Cacciagrano, F. Corradini, E. Merelli, L. Tesei, BioShape: A spatial shape-based scale-independent simulation environment for biological systems, Procedia Computer Science 1 (1) (2010) 827–835 (May 2010).  doi:10.1016/j.procs.2010.04.090.
  * [5] H. V. Westerhoff, M. Verma, F. J. Bruggeman, A. Kolodkin, M. Swat, N. Hayes, M. Nardelli, B. M. Bakker, J. L. Snoep, From Silicon Cellsilicon cell to Silicon Humansilicon human, in: B. Booß-Bavnbek, B. Klösgen, J. Larsen, F. Pociot, E. Renström (Eds.), BetaSys, Springer New York, New York, NY, 2011, pp. 437–458 (2011).  doi:10.1007/978-1-4419-6956-9_19.
  * [6] B. Teusink, J. Passarge, C. A. Reijenga, E. Esgalhado, C. C. van der Weijden, M. Schepper, M. C. Walsh, B. M. Bakker, K. van Dam, H. V. Westerhoff, J. L. Snoep, Can yeast glycolysis be understood in terms of in vitro kinetics of the constituent enzymes? Testing biochemistry: Do we understand yeast glycolysis?, Eur. J. Biochem. 267 (17) (2000) 5313–5329 (Sep. 2000).  doi:10.1046/j.1432-1327.2000.01527.x.
  * [7] N. Cannata, F. Corradini, E. Merelli, L. Tesei, Agent-Based Models of Cellular Systems, in: B. Reisfeld, A. N. Mayeno (Eds.), Computational Toxicology, Vol. 930 of Methods in Molecular Biology, Humana Press, Totowa, NJ, 2013, pp. 399–426 (2013).  doi:10.1007/978-1-62703-059-5_18.
  * [8] E. Merelli, G. Armano, N. Cannata, F. Corradini, M. d’Inverno, A. Doms, P. Lord, A. Martin, L. Milanesi, S. Moller, M. Schroeder, M. Luck, Agents in bioinformatics, computational and systems biology, Briefings in Bioinformatics 8 (1) (2006) 45–59 (May 2006).  doi:10.1093/bib/bbl014.
  * [9] I. Nardecchia, M. Lechelon, M. Gori, I. Donato, J. Preto, E. Floriani, S. Jaeger, S. Mailfert, D. Marguet, P. Ferrier, M. Pettini, Detection of long-range electrostatic interactions between charged molecules by means of fluorescence correlation spectroscopy, Phys. Rev. E 96 (2) (2017) 022403 (Aug. 2017).  doi:10.1103/PhysRevE.96.022403.
  * [10] K. A. Johnson, R. S. Goody, The Original Michaelis Constant: Translation of the 1913 Michaelis–Menten Paper, Biochemistry 50 (39) (2011) 8264–8269 (Oct. 2011).  doi:10.1021/bi201284u.
  * [11] G. E. Briggs, J. B. S. Haldane, A Note on the Kinetics of Enzyme Action, Biochem. J. 19 (2) (1925) 338–339 (Jan. 1925).  doi:10.1042/bj0190338.
  * [12] M. Piangerelli, S. Maestri, E. Merelli, Visualising 2-simplex formation in metabolic reactions, J. Mol. Graph. Model. 97 (2020) 107576 (2020).  doi:10.1016/j.jmgm.2020.107576.
  * [13] H.-J. Kreowski, G. Rozenberg, Graph transformation through graph surfing in reaction systems, Journal of Logical and Algebraic Methods in Programming 109 (2019) 100481 (Dec. 2019).  doi:10.1016/j.jlamp.2019.100481.
  * [14] A. Ehrenfeucht, G. Rozenberg, Reaction systems, Fundam Inf 75 (1–4) (2007) 263–280 (Jan. 2007).
  * [15] S. Maestri, E. Merelli, Algebraic Characterisation of Non-coding RNA, in: P. Cazzaniga, D. Besozzi, I. Merelli, L. Manzoni (Eds.), Computational Intelligence Methods for Bioinformatics and Biostatistics, Vol. 12313, Springer International Publishing, Cham, 2020, pp. 145–158 (2020).  doi:10.1007/978-3-030-63061-4_14.
  * [16] A. Bagossy, G. Vaszil, Simulating reversible computation with reaction systems, J Membr Comput 2 (3) (2020) 179–193 (Oct. 2020).  doi:10.1007/s41965-020-00049-9.
  * [17] A. Ehrenfeucht, I. Petre, G. Rozenberg, Reaction Systems: A Model of Computation Inspired by the Functioning of the Living Cell, WORLD SCIENTIFIC, 2017, pp. 1–32 (Jun. 2017).  doi:10.1142/9789813148208_0001.