Classification of Discrete - Ricerca sull'intelligenza artificiale specializzata al MIT

Classification of Discrete
Dynamical Systems Based
on Transients

In order to develop systems capable of artificial evolution,

Astratto
we need to identify which systems can produce complex behavior.
We present a novel classification method applicable to any class of
deterministic discrete space and time dynamical systems. The method
is based on classifying the asymptotic behavior of the average
computation time in a given system before entering a loop. We were
able to identify a critical region of behavior that corresponds to a phase
transition from ordered behavior to chaos across various classes of
dynamical systems. To show that our approach can be applied to many
different computational systems, we demonstrate the results of
classifying cellular automata, Turing machines, and random Boolean
networks. Further, we use this method to classify 2D cellular automata
to automatically find those with interesting, complex dynamics. Noi
believe that our work can be used to design systems in which complex
structures emerge. Also, it can be used to compare various versions of
existing attempts to model open-ended evolution (Channon, 2006;
Ofria & Wilke, 2004; Ray, 1991).

Barbora Hudcová*
Charles University, Prague
Czech Institute of Informatics,

Robotics, and Cybernetics, Prague
bara.hudcova@gmail.com

Tomáš Mikolov
Czech Institute of Informatics,

Robotics, and Cybernetics, Prague
tmikolov@gmail.com

Keywords
Classification of complex systems,
cellular automata, Turing machines,
random Boolean networks,
transients, phase transition

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

1 introduzione

There are many approaches to searching for systems capable of open-ended evolution. One option
is to carefully design a model and observe its dynamics. Iconic examples were designed by Ray
(1991), Ofria and Wilke (2004), Channon (2006), and Soros and Stanley (2012). Tuttavia, as we
lack any universally accepted formal definition of open-endedness or complexity, there is no formal
method of proving the system is indeed “interesting.” Conversely, lacking definitions of such key
terms, it seems extremely difficult to design such models systematically.

Approaching the problem of searching for open-endedness bottom up, we can define a suitable
classification of dynamical systems that would help us identify a region of complexity. An ideal clas-
sification would be based on a formally defined property, be effectively computable, and help us
automatically search for complex systems possibly capable of modeling artificial evolution.

Over the years, many different metrics have been introduced to study a systemʼs dynamics. As an
esempio, cellular automata (CAs) were studied in terms of their space-time dynamics observations
(Wolfram, 1984), their space-time compression sizes (Zenil, 2010), their actions on probability measures
(Gutowitz, 1990), the Z-parameter (Wuensche & Lesser, 2001), and the lambda parameter (Langton,

* Corresponding author.

Artificial Life 27: 220–245 (2021) https://doi.org/10.1162/artl_a_00342

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

1986). Most of such approaches show that the complex region of systems lies somewhere “in between”
the ordered and chaotic phase.

In questo articolo, we introduce a novel method of classifying complex systems based on estimating
their asymptotic average computation time with increasing space size. The key result is that the clas-
sification identifies a region of systems that seem to be at a phase transition between ordered and
chaotic behavior. Across various classes of discrete systems, we demonstrate that complex systems
such as CAs computing nontrivial tasks, universal Turing machines (TMs), and random Boolean net-
works (RBNs) at a critical phase belong to this region. Even though we are far from characterizing com-
plexity, we hope this method helps us understand which formally defined properties correlate with it.
This article is an extension of Hudcová and Mikolov (2020). Compared to our previous work
where we concentrated solely on CAs, this article includes additional results of the classification
method of TMs and RBNs as well as general conclusions about the average computation time at
a phase transition.

2 Transient Classification: A General Method

We first introduce the basic principle of classification based on transients, which can be applied to
any deterministic discrete space and time dynamical system (DDDS). In subsequent sections, we
describe the results of the classification applied to CAs, TMs, and RBNs to demonstrate its use
across different classes of discrete dynamical systems.

2.1 Basic Notions
Let us consider a generic deterministic discrete system D operating on finite space, characterized by
a tuple D = (S, F ) where S is a finite set of configurations and F : S → S is a global transition
function governing the dynamics of the system. We define the trajectory of a configuration u 2 S as the
sequence
(cid:1)

(cid:3)
tu; F uð Þ; F 2 uð Þ; …

As S is finite, every trajectory eventually becomes periodic. We call the preperiod of this sequence
the transient of initial configuration u and denote its length by tu. More formally, we define tu to be the
smallest positive integer i, for which there exists j 2 N, j > i, such that F i(tu) = F j(tu). The periodic
part of the sequence is called an attractor. The phase-space of D = (S, F) is an oriented graph with
vertices V = S and edges E = {(tu, F(tu)), tu 2 S}. Such a graph is composed of components, each
containing one attractor and multiple transient paths leading to the attractor. The phase-space
completely characterizes the dynamics of the system. Tuttavia, it is infeasible to describe when
the configuration space S is large. Given a DDDS D, we will focus on studying its average transient
length

T Dð Þ ¼

1
Sj j

u2S

tu:

We describe the method of estimating a systemʼs average transient length together with the error
analysis in section 2.4.

2.2 The Main Principle
Suppose we have a sequence of DDDSs

D1 ¼ S1; F1
ð

Þ; D2 ¼ S2; F2
ð

Þ; D3 ¼ S3; F3
ð

Þ; …

Artificial Life Volume 27, Number 3–4

221

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Figura 1. Diagram depicting the asymptotic growth of average transient lengths of a sequence of discrete systems.

operating on configuration spaces of growing size. Questo è, Fi : Si → Si and |Si| < |Si +1| for each i. For instance, the sequence can be given by a cellular automaton with a fixed local rule, operating on a finite cyclic grid of growing size. Our goal is to estimate the asymptotic growth of the systemʼs average transient lengths, as shown in Figure 1. In practice, we generate a finite part of the sequence (|Si|, T(Di)ÞB i¼1 where B is an upper bound imposed by our computational limitations and examine different regression fits of the data. Specifically, we evaluate the fit to constant, logarithmic, linear, polynomial, and exponential functions. We pick the best fit with respect to the R2 score and obtain the classes: Bounded, Log, Lin, Poly, and Exp. If the score of the fit to all such functions is low (i.e., R2 < 85%), we say the system is Unclassified. Surprisingly, we found a very good fit to one of the classes with R2 > 90% for most DDDSs we examined. The trend we
have observed, which seems to hold across various families of DDDSs, is shown in Figure 2. We describe
it in more detail for each family in the subsequent sections.

We do not claim our method determines the true asymptotic behavior of a system; it is merely a
possible interpretation of the method. For some systems, the transient growth might correspond to
more complicated functions but we have deliberately chosen the classes to be quite robust and
coarse to have clearer boundaries between them.

The uncertainty of the true asymptotic growth is especially relevant for the Lin and Poly Classes
which identify the critical phase transition region. Such systems might turn out to be logarithmic or
exponential, and it might be the case that we have not detected this due to our limited data. How-
ever, in such a case, such systems would exhibit significantly slower convergence to their asymptotic
behavior than systems in other classes, which is a typical property of a system at a phase transition.

2.3 Computational Interpretation
In non-classical models of computation (Stepney, 2012), the process of traversing a discrete systemʼs
transients can be perceived as the process of self-organization, in which information can be aggre-
gated in an irreversible manner. The attractors are then viewed as memory storage units, from which

Figura 2. General trend of the transient classification results.

222

Artificial Life Volume 27, Number 3–4

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

the information about the output can be extracted. For CAs this is explored in Kaneko (1986). IL
average transient growth then corresponds to the average computation time of the system1. There-
fore, we can interpret our goal as trying to estimate a systemʼs asymptotic average computation time.
DDDSs with bounded transient lengths can only perform trivial computation. D'altra parte,
DDDSs with exponential transient growth can be interpreted as inefficient computation models.

In the context of artificial evolution, we can view the global transition rule of a DDDS as the phys-
ical rule of the system, whereas the initial configuration is the particular “setting of the universe,” which
is then subject to evolution. If we are interested in finding DDDSs capable of complex behavior
automatically, it would be beneficial for us if such behavior occurred on average, rather than having
to select the initial configurations carefully from some narrow region. The probability of finding such
special initial configurations would be extremely low as the configuration space tends to be very large.
This motivates our study of the growth of average transient lengths rather than the maximum ones.

2.4 Average Transients: Error Estimate
Let us fix a DDDS D = (S, F) operating on a large configuration space, per esempio., |S| ≫ 2100. In such a
case, computing the average transient length (cid:2) is infeasible. Così, we uniformly randomly sample
M
initial configurations u1, u2, …, um and estimate (cid:2) by 1
i¼1 tui. It remains to estimate the number of
M
− (cid:2)| is reasonably small.
samples m so that the error |1
M

More formally, for D = (S, F), let (S, P) be a discrete probability space where S is the set of all
configurations and P is a uniform distribution. Let X : S → N be a random variable, which sends
each u to its transient length tu. This gives rise to a probability distribution of transient lengths on N
with mean E(X) and variance var (X). It can be easily shown that E(X) = (cid:2). Our goal is to obtain a
good estimate of E(X) by the Monte Carlo method (Owen, 2013).

P
M
i¼1 tui

Let (X1, X2, …, Xm) be a random sample of iid random variables, Xi ¼d X for all i. Let (cid:2)(M) = 1

M
i¼1

Xi be the sample mean and

(cid:3) mð Þ ¼

ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Þ2

Xi −(cid:2) mð Þ

i¼1

1
m−1

the sample standard deviation. As X is a mapping from a finite set, var(X) < ∞, and thus we have by the central limit theorem the convergence to a normal distribution. The interval (cid:6) (cid:5) (cid:2) mð Þ − u1 − (cid:4) 2 (cid:3) mð Þ ﬃﬃﬃﬃ p ; (cid:2) mð Þ þ u1 − (cid:4) m 2 (cid:3) mð Þ ﬃﬃﬃﬃ p m where ub is the (cid:5) quantile of the normalized normal distribution covers (cid:2) for m large with probability approximately 1 − (cid:4). We will take (cid:4) = 0.05. Hence, with probability approximately 95% (cid:7) (cid:7) (cid:2) − (cid:2) mð Þ (cid:7) (cid:7) < u0:975 (cid:3) mð Þ ﬃﬃﬃﬃ p : m From the nature of our data, both the values E(X) = (cid:2) and var (X) tend to grow with increasing size of the configuration space. Therefore, to employ a general method of estimating the number of samples, we normalize the error by the sample mean and consider (cid:7) (cid:7) (cid:7) (cid:7) (cid:2)−(cid:2) mð Þ (cid:2) mð Þ : 1 Here, the computation time is understood in the abstract sense; as the number of iterations of the transition function. Artificial Life Volume 27, Number 3–4 223 l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . e d u a r t l / / l a r t i c e - p d f / / / / 2 7 3 – 4 2 2 0 2 0 0 3 3 0 3 a r t l / _ a _ 0 0 3 4 2 p d . f b y g u e s t t o n 0 8 S e p e m b e r 2 0 2 3 B. Hudcová and T. Mikolov Classification of Discrete Dynamical Systems Therefore for m sufficiently large such that u0:975 (cid:3) mð Þ ﬃﬃﬃﬃ p (cid:2) mð Þ m < (cid:6) (1) we have that (cid:2)(m) differs from (cid:2) by at most (cid:6) · 100% with probability approximately 95%. In practice, we put (cid:6) = 0.1 and produce the observations in batches of size 20 until the condition in Equation 1 is met (for most elementary CAs this was satisfied typically after 400 data points were sampled). We approximate the uniform random sampling of initial configurations using Pythonʼs numpy.random library (https://numpy.org/doc/1.16/reference/routines.random.html). 3 Cellular Automata 3.1 Introducing Cellular Automata Informally, a cellular automaton can be perceived as a k-dimensional grid consisting of identical finite state automata. They are all updated synchronously in discrete time steps based on an identical local update function depending only on the states of automata in their local neighborhood. A for- mal definition can be found in Jarkko Kari (2005). CAs were first studied as models of self-replicating structures (Langton, 1984; Neumann & Burks, 1966; Reggia et al., 1993). Subsequently, they were examined as dynamical systems (Hedlund, 1969; Gutowitz et al., 1987; Vichniac, 1984), or as models of computation (Mitchell, 1998; Toffoli, 1977). Being so simple to simulate, yet capable of complex behavior and emergent phenomena (Crutchfield & Hanson, 1993; Hanson, 2009), CAs provide a convenient tool to examine the key, yet undefined notions of complexity and emergence. 3.1.1 Basic Notions We study the simple class of elementary cellular automata (ECAs), which are one-dimensional CAs with two states {0, 1} and neighborhood of size 3. Each ECA is given by a local transition function f :{0, 1}3 → {0, 1}. Hence, there are only 256 of them. The size of this CA class is the reason for making it our first case study. One can simply explore it by studying the dynamics of every single ECA. We identify each local rule f determining an ECA with Wolfram number of the ECA defined as: 20f 0; 0; 0 ð Þ þ 21f 0; 0; 1 ð Þ þ 22f 0; 1; 0 ð Þ þ … þ 27f 1; 1; 1 ð Þ: We will refer to each ECA as a “rule k” where k is the corresponding Wolfram number of its un- derlying local rule. We will consider the CA to operate on finite grids with periodic boundary conditions. Hence, given a local rule f and a grid size n, we obtain a configuration space {0, 1}n and a global update function F : {0, 1}n → {0, 1}n. Let ({0, 1}n, F) be an ECA operating on a grid of size n and (u, F(u), F2(u), …) a trajectory of a configuration u 2 {0, 1}n. The space-time diagram of such a simulation is obtained by plotting the configurations as horizontal rows of black and white squares (corresponding to states 1 and 0) with a vertical axis determining the time, which is progressing downwards. We note that properties of CA phase-spaces were examined by, among others, Wuensche and Lesser (2001). Precisely for this purpose, a software was designed by Wuensche (2016). 224 Artificial Life Volume 27, Number 3–4 l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . e d u a r t l / / l a r t i c e - p d f / / / / 2 7 3 – 4 2 2 0 2 0 0 3 3 0 3 a r t l / _ a _ 0 0 3 4 2 p d . f b y g u e s t t o n 0 8 S e p e m b e r 2 0 2 3 B. Hudcová and T. Mikolov Classification of Discrete Dynamical Systems 3.2 History of CA Classifications An ideal classification would be based on a rigorously defined and efficiently measurable property, identifying a region of systems with interesting behavior. In this section, we describe three qualita- tively different classifications of ECAs, and subsequently, we will compare our results to them. 3.2.1 Wolframʼs Classification The most intuitive and simple approach to examining the dynamics of CAs is to observe their space- time diagrams. This method was particularly proclaimed by Wolfram (2002). Therein, he established an informal classification of CA dynamics based on such diagrams. Wolfram (2002, p. 231–249) distinguishes the following classes, which are shown in Figure 3. Class 1 … quickly resolves to a homogenous state Class 2 … exhibits simple periodic behavior Class 3 … exhibits chaotic or random behavior Class 4 … produces localized structures that interact with each other in complicated ways The main issue is that we have no formal method of classifying CAs in this way. In fact, this problem is in general undecidable (Culik & Yu, 1988). Moreover, the behavior of some CAs can vary with l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . e d u a r t l / / l a r t i c e - p d f / / / / 2 7 3 – 4 2 2 0 2 0 0 3 3 0 3 a r t l / _ a _ 0 0 3 4 2 p d . f b y g u e s t t o n 0 8 S e p e m b e r 2 0 2 3 Figure 3. Space-time diagrams of rules from each of Wolframʼs (2002) classes. Class 1 rule 32 is on top left, Class 2 rule 108 on top right, Class 4 rule 110 on the bottom left, and Class 3 rule 30 on the bottom right. Artificial Life Volume 27, Number 3–4 225 B. Hudcová and T. Mikolov Classification of Discrete Dynamical Systems Figure 4. On the left, rule 126 is simulated with an initial condition consisting of a single 1 bit padded with 0ʼs. On the right, the same rule is simulated with a random initial configuration. different initial configurations. An example being rule 126, which oscillates between Class 2 and Class 3 behavior, as shown in Figure 4. The transient classification we present in this article deals with both these issues. 3.2.2 Zenilʼs Classification In the first part of his article Zenil (2010) studied the compression size of the space-time diagrams of each ECA simulated for a fixed amount of steps. For the classification, he examines the simu- lations from a particular initial configuration (a single one surrounded by zeros). Using a clustering technique, Zenil obtained two classes distinguishing between Wolframʼs simple classes 1 and 2 and com- plex classes 3 and 4. We show our reproduction of Zenilʼs results in Figure 5. Zenilʼs method nicely formalizes Wolframʼs observations of the space-time diagrams. However, the results depend on the choice of initial conditions as well as the grid size, data representation, and the compression algorithm. We conducted multiple experiments presented in Figure 6, which suggest that Zenilʼs results might be sensitive to the choice of such parameters. We note that he addresses the sensitivity to the choice of the initial configurations in the second part of his article (Zenil, 2010). In vast CA spaces where it is not feasible to examine every CA and mark it into one of Wolframʼs classes by hand, it is not clear how the parameter values should be chosen. Moreover, the data rep- resentation causes the extension of this method to more general dynamical systems to be problematic; for example, using gzip (https://www.gnu.org/software/gzip/) to compress space-time diagrams of a 2D CA is suboptimal. Figure 5. Reproduction of Zenilʼs results (Zenil, 2010). The purple cluster corresponds to the interesting Class 3 and 4 rules, the yellow cluster to the rest. 226 Artificial Life Volume 27, Number 3–4 l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . e d u a r t l / / l a r t i c e - p d f / / / / 2 7 3 – 4 2 2 0 2 0 0 3 3 0 3 a r t l / _ a _ 0 0 3 4 2 p d . f b y g u e s t t o n 0 8 S e p e m b e r 2 0 2 3 B. Hudcová and T. Mikolov Classification of Discrete Dynamical Systems Figure 6. Graphs representing the results of Zenilʼs (2010) method when different parameter values are used that dem- onstrate the possible sensitivity of the results. On the left, the ECAs were simulated for a longer time, which caused complex rules 110, 124, 137, and 193 no longer to belong to the interesting purple cluster. On the right, the ECAs were simulated from a fixed, randomly chosen initial condition. In such a case, we obtain entirely different clusters. ECA = elementary cellular automata. 3.2.3 Wuenscheʼs Z-parameter Wuensche and Lesser (2001) chose an interesting approach by studying ECA behavior when reversing the simulations and computing the preimages of each configuration. They introduce the Z-parameter, representing the probability that a partial preimage can be uniquely prolonged by one symbol and sug- gest that Class 4 CAs typically occurs at Z ≈ 0.75. However, no clear classification is formed. The crucial advantage is that the Z-parameter depends only on the CAʼs local rule and can be com- puted effectively. It is, however, questionable whether studying only the local rule could describe the overall dynamics of a system sufficiently well. We note that transients of CAs have been examined, as in Wuensche and Lesser (2001) and Gutowitz (1994). However, we are not aware of an attempt to compare the asymptotic growth of transients for different ECAs. 3.3 Transient Classification of ECAs For each ECA given by a local rule f, we consider the sequence of systems (cid:3) (cid:3) (cid:1) D3 ¼ 0; 1f g3; F3 (cid:1) ; D4 ¼ 0; 1f g4; F4 ; … l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . e d u a r t l / / l a r t i c e - p d f / / / / 2 7 3 – 4 2 2 0 2 0 0 3 3 0 3 a r t l / _ a _ 0 0 3 4 2 p d . f b y g u e s t t o n 0 8 S e p e m b e r 2 0 2 3 Figure 7. Bounded Class rule 36. The average transient plot is on the left, the space-time diagram on the right. Artificial Life Volume 27, Number 3–4 227 B. Hudcová and T. Mikolov Classification of Discrete Dynamical Systems Figure 8. Log Class rule 28. The average transient plot is on the left, the space-time diagram on the right. that represents the ECAs operating on grids of growing size. We can apply the transient classifica- tion to this sequence, as described in section 2, to estimate the asymptotic growth of the average transient lengths for each ECA. We consider all 256 ECAs up to equivalence classes obtained by changing the role of left and right neighbor, the role of 0 and 1 state, or both. It can be easily shown that automata in the same equivalence class have isomorphic phase spaces for any grid size. Thus, they perform the same com- putation. This yields 88 effectively different ECAs, each being a representative with the minimum Wolfram number from its corresponding equivalence class. In this section, we present the classifi- cation of the 88 unique ECAs based on their asymptotic transient growth. 3.3.1 Results We obtained a surprisingly clear classification of all 88 unique ECAs with four major classes corre- sponding to the bounded, logarithmic, linear, and exponential growth of average transients. Below, we give a more detailed description of each class. Bounded Class: 27/88 rules (30.68%) The average transient lengths were bounded by a constant independent of the grid size. This suggests that the long-term dynamics of such automata can be predicted efficiently. See Figure 7. Figure 9. Lin Class rule 62. The average transient plot is on the left, the space-time diagram on the right. 228 Artificial Life Volume 27, Number 3–4 l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . e d u a r t l / / l a r t i c e - p d f / / / / 2 7 3 – 4 2 2 0 2 0 0 3 3 0 3 a r t l / _ a _ 0 0 3 4 2 p d . f b y g u e s t t o n 0 8 S e p e m b e r 2 0 2 3 B. Hudcová and T. Mikolov Classification of Discrete Dynamical Systems Figure 10. Exp Class rule 45. The average transient plot is on the left, the space-time diagram on the right. Log Class: 39/88 rules (44.32%) The largest ECA class exhibits logarithmic average transient growth. The event of two cells at a large distance “communicating” is improbable for this class. See Figure 8. Lin Class: 8/88 rules (9.09%) On average, information can be aggregated from cells at an arbitrary distance. This class contains automata whose space-time diagrams resemble some sort of computation. This is supported by the fact that this class contains two rules known to have a nontrivial computa- tional capacity: rule 184, which computes the majority of black and white cells, and rule 110, which is the only ECA so far proven to be Turing complete (Cook, 2004). We note that rules in this class are not necessarily complex as the interesting behavior seems to correlate with the slope of the linear growth. Most of Lin Class rules had only a very gradual incline. In fact, the only two rules with such slope greater than 1, rules 110 and 62, seem to be the ones with the most interesting space-time diagrams. See Figure 9. We are aware that average transients of rules in Lin Class might turn out to grow logarithmically or exponentially given enough data samples. In such a case, the rules in Lin Class show a significantly slower convergence to their asymptotic behavior, which supports the hypothesis that they belong to a phase transition region. Exp Class: 6/88 rules (6.82%) This class has a striking correspondence to automata with chaotic behavior. Visually, there seem to be no persistent patterns in the configurations. Not only the tran- sients but also the attractor lengths are significantly larger than for other rules. The rules with the fastest growing transients are 45, 30, and 106. See Figure 10. Figure 11. Affine Class rule 90. The average transient plot is on the left, the space-time diagram on the right. Artificial Life Volume 27, Number 3–4 229 l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . e d u a r t l / / l a r t i c e - p d f / / / / 2 7 3 – 4 2 2 0 2 0 0 3 3 0 3 a r t l / _ a _ 0 0 3 4 2 p d . f b y g u e s t t o n 0 8 S e p e m b e r 2 0 2 3 B. Hudcová and T. Mikolov Classification of Discrete Dynamical Systems Figure 12. Fractal Class rule 126. The average transient plot is on the left, the space-time diagram on the right. Affine Class: 4/88 rules (4.55%) This class contains rules 60, 90, 105, and 150 whose local rules are affine Boolean functions. Such automata can be studied algebraically and predicted efficiently. It was shown in Martin et al. (1984) that the transient lengths of rule 90 depend on the largest power of 2, which divides the grid size. Therefore, the measured data did not fit any of the functions above but formed a rather specific pattern. See Figure 11. Fractal Class: 4/88 rules (4.55%) This class contains rules 18, 122, 126, and 146 which are sensitive to initial conditions. Their evolution either produces a fractal structure resembling a Sierpinski triangle or a space-time diagram with no apparent structures. We could say such rules oscillate between easily predictable behavior and chaotic behavior. Their average transients and periods grow quite fast, making it difficult to gather data for larger grid sizes. See Figure 12. 3.4 Discussion We have also tried to measure the asymptotic growth of the average attractor size au, u 2 {0, 1}n as well as the average rho value defined as (cid:7)u = tu + au. This, however, produced data points that could not be fitted to simple functions well. This is due to the fact that many automata have attractors consisting of a configuration, which is shifted by one bit to the left, resp. right, at every time step. The size of such an attractor then depends on the greatest common divisor of the size of the period of the attractor and the grid size, and this causes oscillations. We conclude that such phase-space properties are not suitable for this classification method. Below, we compare our results to other classifications described earlier. An exhaustive compar- ison for each ECA is presented in Figure 13. Wolframʼs Classification. The significance of our results for ECA stems precisely from the fact that the transient classification corresponds to Wolframʼs (2002) so well. As it is not clear for many rules which Wolfram class they belong to, the main advantage is that we provide a formal criterion upon which this could be decided. In particular, rules in Bounded Class and Log Class correspond to rules in either Class 1 or Class 2. Exp Class corresponds to the chaotic Class 3, and Lin Class contains Class 4 together with some Class 2 rules. We note an interesting discrepancy: Rule 54, which is possibly considered by Wolfram to be Turing complete, belongs to the Exp Class. This might suggest that computations performed by this rule can be on average quite inefficient. Zenilʼs Classification. Zenilʼs (2010) classification of ECAs offers a great formalization of Wolframʼs and seems to roughly correspond to it. Compared to the transient classification, it is, however, less fine-grained. Moreover, it contains some arbitrary parameters, such as the data representation and compression algorithm used. In addition, it uses a clustering technique, which requires data of multiple 230 Artificial Life Volume 27, Number 3–4 l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . e d u a r t l / / l a r t i c e - p d f / / / / 2 7 3 – 4 2 2 0 2 0 0 3 3 0 3 a r t l / _ a _ 0 0 3 4 2 p d . f b y g u e s t t o n 0 8 S e p e m b e r 2 0 2 3 B. Hudcová and T. Mikolov Classification of Discrete Dynamical Systems l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . e d u a r t l / / l a r t i c e - p d f / / / / 2 7 3 – 4 2 2 0 2 0 0 3 3 0 3 a r t l / _ a _ 0 0 3 4 2 p d . f b y g u e s t t o n 0 8 S e p e m b e r 2 0 2 3 Figure 13. Comparing classifications of the 88 unique ECAs. Adapted from Wolfram (2002), Zenil (2010), and Wuensche and Lesser (2001). automata to be mutually compared in order to give rise to different classes. In contrast, the transient class can be determined for a single automaton without any context. Another important difference is that Zenil observed the simulations from a fixed initial config- uration; that is, he examined the local dynamics of ECAs. In contrast, the transient classification studies their global dynamics. Wuenscheʼs Z-parameter. Wuensche (Wuensche & Lesser, 2001) suggests that complex behavior occurs around Z = 0.75, which agrees with the fact that Lin Class rules with a steep slope (rule 110, 62, and 25) have this Z value precisely. However, the Z = 0.75 is in fact quite frequent. This sug- gests that thanks to its simplicity, the Z parameter can be used to narrow down a vast space of CA rules when searching for complexity. However, more refined methods have to be subsequently ap- plied to find concrete CAs with interesting behavior. Artificial Life Volume 27, Number 3–4 231 B. Hudcová and T. Mikolov Classification of Discrete Dynamical Systems Table 1. Classification of 10,000 randomly sampled symmetric 2D 3-state CAs. Transient Class Bounded Log Lin Poly Exp Unclassified Percentage of CAs 0% 18.21% 1.17% 1.03% 72.62% 6.97% 3.5 Transient Classification of 2D CAs So far, we have examined the toy model of ECAs. Transient classificationʼs true usefulness would stem from its application to more complex CAs, where it could be used to discover automata with interesting behavior. Therefore, we applied the classification on a subset of two-dimensional CAs with a 3 × 3 neighborhood and three states to see whether 2D automata would still exhibit such clear transient growths. We work with 2D CAs operating on a finite square grid of size n × n. We consider the topology of the grid to be that of a torus for each cell to have a uniform neighborhood. We estimated the average transient length and measured the asymptotic growth with respect to n (i.e., the size of the square gridʼs side). This is motivated by the fact that in an n × n grid, the greatest distance between two cells depends linearly on n rather than quadratically. To reduce the vast automaton space, we only considered such automata whose local rules are invariant to all the symmetries of a square. As there are still 32861 such symmetrical 2D CAs, we randomly sampled 10,000 of them. For such a large space, we could not examine each CA individually. Therefore, we fit the average transient growth to bounded, logarithmic, linear, polynomial, and exponential functions to obtain the Bounded, Log, Lin, Poly, and Exp classes. If none of the fits gave a good enough score (i.e., R2 >
85%), then we marked the corresponding CA as unclassified. We were able to classify 93.03% Di
10,000 sampled automata with a time bound of 40 s for the computation of one transient length
value on a single CPU. We estimate that most CAs are unclassified due to such a computation re-
sources restriction, or rather strict conditions we imposed on a good regression fit. In this large
space of 2D CAs Exp Class seems to dominate the rule space. Another interesting aspect in which
2D CAs differ from the ECAs is the emergence of rules in Poly Class; the transients of such rules
grow approximately quadratically. Inoltre, our results suggest that the occurrence of Bounded
Class CAs in 2D is much scarcer as we found no such CAs in our sample. See Table 1.

We observed the space-time diagrams of randomly sampled automata from each class to infer its
typical behavior. On average, Log Class automata quickly enter attractors of small size. Lin Class
exhibits the emergence of various local structures. For automata with a more gradual incline, come
structures seem to die out quite fast. Automata with steeper slopes exhibit complex interactions of
such structures. Poly Class automata with a steep slope seem to produce spatially separated regions
of chaotic behavior against a static background. In the case of more gradual slopes, some local struc-
tures emerge. Finalmente, Exp Class seems to be evolving chaotically with no apparent local structures.
We present various examples of CA evolution dynamics in the form of GIF animations here2.

This suggests that the regions of Lin Class with a steep slope and Poly Class with a more gradual
incline seems to contain a non-trivial ratio of automata with complex behavior. In this sense, IL

https://bit.ly/trans_class

232

Artificial Life Volume 27, Number 3–4

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Figura 14. Game of Life (Gardener, 1970). The average transient growth plot is on the left. On the right, we show a
space-time diagram at time t = 200 started from a random initial configuration.

transient classification can assist us to automatically search for complex automata similarly to the
method designed by Cisneros et al. (2019) where interesting novel automata were discovered by
measuring growth of structured complexity using a data compression approach.

3.6 Transients Classification of Other Well-Known CAs
We were interested in whether some well-known complex automata from larger CA spaces would
conform to the transient classification as well. As we show in this section, the result is positive.
Game of Life. As the left plot in Figure 14 suggests, the Turing complete Game of Life (Gardener,
1970) seems to fit Lin Class. This is confirmed by the linear regression fit with R2 ≈ 98.4%.
Genetically Evolved Majority CAs. Mitchell et al. (2000) studied how genetic algorithms can
evolve CAs capable of global coordination. The authors were able to find a 1D CA denoted as (cid:8)par
with two states and radius r = 3, which is quite successful at computing the majority task with the
output required to be of the form of a homogenous state of either all 0ʼs or all 1ʼs. This CA seems to
belong to Lin Class, which is confirmed by the linear regression fit with R2 ≈ 99.2%. See Figure 15.
Totalistic 1D 3-state CAs. A totalistic CA is any CA whose local rule depends only on the number
of cells in each state and not on its particular position. Wolfram studied various CA classes, one of
them being the totalistic 1D CA with radius r = 1 E 3 states S = {0, 1, 2}.

Figura 15. Cellular automaton (CA) (cid:8)par. The average transient growth plot is on the left. On the right, we show a space-
time diagram simulated from a random initial configuration.

Artificial Life Volume 27, Number 3–4

233

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Figura 16. Totalistic cellular automaton with code 1635. The average transient growth plot is on the left. On the right,
we show a space-time diagram of the evolution from a random initial configuration.

Wolfram (2002) presents a list of possibly complex CAs from this class. We applied the transient
classification to such CAs and learned that most of them were classified as logarithmic. This agrees
with our space-time diagram observations that the local structures in such CAs “die out” quite quickly.
Nonetheless, some of the CAs were classified as linear. An example of such a CA is in Figure 16 Dove
the linear regression fit has R2 ≈ 97.63%.

4 Turing Machines

In order to demonstrate the generality of the transient classification method, we further used it to
examine the dynamics of TMs. In this section, we present the classification results.

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

4.1 Introducing Turing Machines
Informally, a TM consists of an infinite tape divided into cells and a movable reading head that scans
one cell of the tape at a time. Every cell contains a symbol from some finite alphabet A, and the TM
is at an internal state from a finite set S. Depending on the symbol the head is reading and on its
internal state, the TM changes its internal state, rewrites the symbol on the tape, and either moves
one cell to the left or right, or stays in place. TMs represent the most classical model of computation;
the Church–Turing thesis states that “effectively calculable functions” are exactly those that can be
realized by a TM (Turing, 1937). For a formal definition of TMs as well as a great introduction to
computability theory, see Soare (2016).

In questo articolo, we will consider deterministic TMs with one tape. S will always denote the finite set
of internal states, A will denote the finite set of tape symbols. To ensure the TM operates on a finite
grid, as in the case of CAs, we will consider a tape of finite size with periodic boundary conditions.
Therefore, each TM with S and A operating on a tape of size n gives rise to a global update function

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

F : An (cid:2) 1; 2; …; N
F

G (cid:2) S → An (cid:2) 1; 2; …; N

G (cid:2) S;

where each configuration specifies the content of the tape, the position of the head, and the internal
state. Così, we can apply the transient classification to it. We emphasize the non-traditional notion
of the halting computation that we consider here. Classically, a TM is considered to halt when it
enters an attractor of size 1; questo è, when it does not change the tapeʼs content, the headʼs position,

234

Artificial Life Volume 27, Number 3–4

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Figura 17. Space-time diagrams of a 6 symbol, 5 state Turing machine in Exp Class. Classical space-time diagram is shown
on the left. On the right we show the space-time diagram depicting the content of the tape after every 50 steps of
computation on a tape of size 50.

or its internal state anymore. In our case, using the interpretation

transients ≈ computation
attractors ≈ memory

we consider a TM to halt whenever it enters any attractor. This is a much weaker notion of halting.
We will depict the space-time diagrams of TM computation as a matrix, each row corresponds to
the content of the tape at subsequent time steps, and time is progressing downwards. As opposed to
CAs with their inherently parallel nature, TMs are sequential computational models. Così, at each
time step, only one symbol on the tape is changed. To produce space-time diagrams comparable to
those produced by CAs, we only depict tape contents at every nth step where n is the size of the tape.
This helps us to intuitively recognize the chaotic dynamics of TMs; Guarda la figura 17. To the best of our
knowledge, we know of no prior work examining the transients of TMs operating on cyclic tapes.

4.2 Transient Classification of Turing Machines
We have studied “small” TMs with the number of states |S| ranging from 4 A 8 and the number of
alphabet symbols |UN| ranging from 2 A 5. For every such combination of values |S| E |UN|, we
have randomly generated 100 transition functions of TMs and computed each TMʼs average tran-
sient length estimate for cyclic tapes of sizes ranging from 20 A 400.

Tavolo 2. Classification of 100 randomly sampled TMs with 7 states and 4
symbols.

Transient Class

Bounded

Log

Lin

Poly

Exp

Unclassified

Percentage of TMs

41%

28%

13%

15%

Artificial Life Volume 27, Number 3–4

235

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Figura 18. Example of a Turing machine with 7 stati, 4 symbols in Lin Class. Its space-time diagram seems to exhibit
nontrivial behavior.

4.2.1 Results
For all the considered values of |S| E |UN|, more than 90% of the TMs were successfully clas-
sified using the transient classification method. We discuss a particular example in more detail below.
Touring Machines with 7 states and 4 symbols. As an example, we present classification results of
100 TMs with 7 states and 4 tape symbols. The results are summarized in Table 2.
Bounded Class: (41/100 TMs) In this space of rather small TMs, Bounded Class seems to dominate
the space. TMs in Bounded Class halt in time independent of the tape size. Therefore, for such TMs it
seems improbable to perform any nontrivial computation on both the finite and infinite tape.
Log Class: (2/100 TMs) Log Class seems to be relatively small across all the TM classes we have
examined. Here, the event that a TM head will read the whole input from the tape is improbable for
large tape sizes.
Lin Class: (28/100 TMs) Let us consider a TM with trivial dynamics, Quale, given any input
configuration, traverses each cell one by one and changes the state of each cell to the state 0 2 S.
After all cells enter this state, the computation halts. Such trivial behavior could be realized in constant
time by a CA, but for a TM it takes at least n steps where n is the size of the tape. Hence, some TMs in
Lin Class exhibit periodic or simple dynamics. This emphasizes the fact that being contained in Lin
Class or Poly Class seems to be a necessary condition for complexity, not a sufficient one. Vedere

Figura 19. Example of a Turing machine with 7 stati, 4 symbols in Poly Class.

236

Artificial Life Volume 27, Number 3–4

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Figura 20. Example of a Turing machine with 7 stati, 4 symbols in Exp Class.

Figura 18. Nevertheless, we have observed TM in the Lin Class whose space-time diagrams seem to
contain some higher-level structures.
Poly Class: (13/100 TMs) In Poly Class we have also observed TMs producing some higher-order struc-
tures. See Figure 19.
Exp Class: (15/100 TMs) We find it interesting that once only every nth row of the space-time
diagram is depicted (n being the tape size). The space-time diagrams of TMs in Exp Class resemble
the space-time diagrams of chaotic CAs. See Figure 20.

As in the case of CAs we are aware that the true asymptotic behavior of TMs in Lin Class or Poly
Class might turn out to be logarithmic or exponential. In such a case, the systems in these classes
would need significantly longer time to converge to their typical long-term behavior, which is a typical
property of systems at a phase-transition.

4.3 Transient Classification of Universal TM
Without much doubt, universal TMs are considered complex. We have estimated the asymptotic
average computation time of 7 universal TMs with a small number of states and symbols constructed
by Rogozhin (1996). All were successfully classified, 6 belonging to Poly Class and 1 to Lin Class. An
example is shown in Figure 21.

Such results agree with the ones obtained for CAs and support the hypothesis that complex

dynamical systems belong to Lin Class or Poly Class.

Figura 21. Universal Turing machines with 10 states and 3 symbols belonging to Poly Class. Again, every 50th step of the
computation is shown in the space-time diagram.

Artificial Life Volume 27, Number 3–4

237

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

5 Random Boolean Networks

Random Boolean networks form a very wide class of discrete dynamical systems that contains both
CAs and TMs. In this section, we show that dynamical systems from this general class also con-
formed to our classification method.

5.1 Introducing Random Boolean Networks
The classical N−K RBN is given by an oriented graph with N nodes, each having exactly K edges
pointing toward it. Inoltre, each node is equipped with a Boolean function of K variables. Every
node can have the value of either 0 O 1; Perciò, the configuration space is exactly {0, 1}N. To update
a particular configuration of the network, the values of all nodes are changed in parallel, according to the
outputs of their corresponding Boolean functions. This gives rise to a global update rule F : {0, 1}N →
{0, 1}N. For a concise introduction to RBNs see Gershenson (2004).

RBNs were first introduced by Stuart Kauffman (1969) as models of gene regulatory networks.
Classically, the nodes are interpreted as genes of a particular organism; their value represents whether
the gene is “turned on” or “off.” In this setting, the attractors of the network represent different cell
types of the organism. Over the years, the networks have been widely studied as models of cell
differentiation (Huang & Ingber, 2000), immune response (Kauffman & Weinberger, 1989), E
neural networks (Wuensche, 1996). The measures of criticality in RBNs have been studied by
Luque and Solé (2000), Wang et al. (2011), and Pineda et al. (2019). For a great overview of work
done on RBNs see Aldana et al. (2003).

5.1.1 Critical Behavior in RBNs
RBNs are generic in the sense that both the connections of nodes and the Boolean functions are
chosen uniformly at random. This makes it possible to analytically study the properties of a typical
N−K network. Infatti, different approaches (Derrida & Pomeau, 1986; Luque & Solé, 1997) Guida
to the same description of phase transitions in RBNs. We describe the results briefly below, as we
will use them in our experiments.

We will describe a slightly more general model of RBNs. We consider a non-uniform connectivity
of the nodes—for a network of size N, we will assign to each node i 2 {1, …, N} the connectivity Ki
2 N and a Boolean function fi of arity Ki. Such a network is parametrized by the mean connectivity
N
hKi = 1
i¼1 Ki. We also introduce the Boolean function sampling bias p 2 (0, 1). Questo è, we will
N
sample the Boolean functions so that for all i the probability that fi(x1, …, xKi ) = 1 is p. Derrida and
Pomeau (1986) have analytically determined the edge of chaos for RBNs depending on the mean
connectivity parameter hKi and Boolean function bias p. By studying the evolution of the distance
between two randomly generated initial configurations over time, they have shown that the critical
values of hKi and p are exactly those satisfying

Kh i ¼

1
2P 1 − p
ð

They obtain the following phases of RBN behavior.

(2)

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

Ordered Phase … RBNs with Kh i < Critical Phase … RBNs with Kh i ¼ Chaotic Phase … RBNs with Kh i >

1
2P 1 − p
ð
1
2P 1 − p
ð
1
2P 1 − p
ð

The curve given by (1) is shown in Figure 22.

238

Artificial Life Volume 27, Number 3–4

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Figura 22. Red curve depicts the critical values of random Boolean network (RBN) mean connectivity hKi and Boolean
function bias p. The blue area denotes the region of ordered behavior; the white area denotes the chaotic region.

In section 5.2, we support the analytical results by showing that the transient classification clearly

distinguished between the ordered, critical, and chaotic regions.

5.1.2 Phase-Space Properties of RBNs
The generic nature of N−K RBNs makes it possible to analytically study their global dynamics. Questo
is a key difference between CAs and RBNs: CAs have a very particular architecture with only local
connections and a uniform local transition rule. Therefore, the mean-field approximation methods
of phase-space properties used for a “typical” N−K RBN would not be as easy to apply to CAs.
With the classical interpretation of RBNs as gene regulatory networks where attractors represent
different cell types, most of the focus has been on analyzing the number and size of the attractors
for different values of N and K. We briefly summarize some results related to our experiments below.
For a more detailed discussion see Aldana et al. (2003).
Ordered Phase. In the ordered phase, when K = 1, it has been shown that a probability of having
an attractor of size l falls exponentially with l (Flyvbjerg & Kjaer, 1988). For a subset of K = 2 RBNs
with ordered behavior, Lynch (1986) has shown that their average transient time grows at most
logarithmically with the network size N.
2 and K ≥ N, the RBN is essentially a random mapping whose
Chaotic Phase. In the case when p = 1
phase-space properties have been studied extensively. It has been shown that both the average
attractor and transient lengths of such RBNs grow exponentially with increasing N (Harris, 1960;
Derrida & Flyvbjerg, 1987).

We note that some previous work examining the transients of RBNs was conducted by Bhattacharjya
and Liang (1996) and Wuensche (1996) but we are not aware of any studies that would use the
asymptotic transient growth to describe the behavior of RBNs at the critical region.

5.2 Transient Classification of RBNs
We have sampled RBNs parametrized by the mean connectivity hKi and the Boolean function bias p.
In this section we show that the results of transient classification clearly distinguish the ordered,
critical, and chaotic phase of RBNs.

5.2.1 Details of the Experiment
Our goal is to estimate the average transient length of a “typical” RBN of size N, with mean connec-
tivity hKi and Boolean function bias p. We would do so for increasing N to observe the asymptotic
behavior. We proceed as follows:

Artificial Life Volume 27, Number 3–4

239

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Figura 23. Growth of typical average transient lengths for RBNs in the ordered region. RBNs with mean connectivity hKi =
1 and Boolean function bias p = 0.5 on the left RBNs with hKi = 4.556 and p = 0.9 on the right. The best fit for both was
logarithmic, with R2 score over 90%.

1. Given N, hKi, and p, we generate a RBN R(N, hKi, P) with the corresponding parameters.
We estimate the average transient length T(R(N, hKi, P)) of R(N, hKi, P) using the
approach described in section 2.4.

2. We repeat step 1 and generate a sequence of RBNs

R1 N ; Kh i; P

Þ; R2 N ; Kh i; P

Þ; …; Rm N ; Kh i; P
ð

and their average transient lengths

T R1 N ; Kh i; P

Þ; …; T Rm N ; Kh i; P
ð

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

to ensure that we are close to the average transient length of an average RBN with
parameters N, hKi, and p. We determine the number of sampled RBNs needed to get
sufficiently close to the true average behavior by a method analogous to the one described

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

Figura 24. Growth of typical average transient lengths for RBNs in the critical region. RBNs with mean connectivity hKi =
2 and Boolean function bias p = 0.5 on the left; RBNs with hKi = 2.381 and p = 0.7 on the right. The best fit for both was
linear, with R2 score over 95%.

240

Artificial Life Volume 27, Number 3–4

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Figura 25. Growth of typical average transient lengths for RBNs in the chaotic region. RBNs with mean connectivity K =
4 and Boolean function bias p = 0.5 on the left; RBNs with K = 4.083 and p = 0.4 on the right. The best fit for both was
exponential, with R2 score over 99%.

in section 2.4. Finalmente, we obtain the typical average transient length as

T N ; Kh i; P

Þ ¼

1
M

T Ri N ; Kh i; P

Þ:

i¼1

3. We try to approximate the sequence (T(N, hKi, P)Þ

typically compute (T(N, hKi, P)Þ200
imposed by the computation time of the transient lengths.

∞
N ¼1 by generating a finite part of it. Noi
N ¼5, the upper bound being either N = 200 or the limit

5.3. Results
Ordered Phase. We have computed Kc , p along the curve given by (1) for p = 0.1, 0.2, …, 0.9 E
sampled RBNs with parameters hKc – 1i, p to ensure we are in the ordered region. See Figure 23.
For all such ensembles of RBNs the best fit for the typical average transient asymptotic growth

was logarithmic. This supports the analytical results proven for special cases of K and p values.
Critical Phase. We have sampled RBNs with parameters hKci, p along the curve given by (1) for p =
0.1, 0.2, 0.3, …, 0.9. In all the sampled cases, the best fit for the typical average transient growth was
linear. As in the case for CAs, we are aware that the asymptotic behavior of such RBNs can turn out
to be logarithmic or exponential and that we just might not have sampled large enough networks. In
such a case, we can interpret Lin Class as a region of RBNs that take significantly longer to converge
to their asymptotic behavior. See Figure 24.
Chaotic Phase. We have computed the critical values Kc , p along the curve given by (1) for values
p = 0.2, 0.3, …, 0.7, 0.8 and sampled RBNs with parameters hKc + 2io, p to ensure we are in the
chaotic region. For all such ensembles of RBNs, the best fit for the typical average transient asymp-
totic growth was exponential, which again agrees with the analytic results. See Figure 25.

These experiments support the results obtained for CAs and TMs indicating that ordered discrete
systems belong to Bounded Class and Log Class, chaotic systems correspond to Exp Class, E
complex systems lie in the region in between, corresponding to Lin Class and Poly Class.

6 Conclusione

We presented a classification method based on the asymptotic growth of average computation time.
It is applicable to any deterministic discrete space and time dynamical system. We presented a good

Artificial Life Volume 27, Number 3–4

241

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

correspondence between transient classification and Wolframʼs (2002) classification in the case of
ECAs. Further, we showed that the classification works for 2D CAs, TMs, and RBNs, and we used
it to discover 2D CAs capable of emergent phenomena. By demonstrating that complex discrete
systems (per esempio., Game of Life, rule 110, and several universal TMs and RBNs with critical parameter
values) belong to Lin Class or Poly Class, we have reason to believe that linear and polynomial tran-
sient growth navigates us toward a region of complex discrete systems.

Another elegant alternative would be to merge Bounded Class and Log Class representing the
ordered phase, and Lin Class and Poly Class corresponding to the critical phase. In this way, we
would obtain the traditional three phases of dynamics. This is entirely possible; we have kept the
five classes to respect our initial experiments on elementary CAs where we obtained a much finer
classification scheme.

The classification is based on a very simple idea and can be implemented with a few lines of code.
In the field of ALife where novel discrete dynamical systems are designed as possible models of
artificial evolution, the method we presented can be used to check whether such systems belong
to Lin Class or Poly Class. This might support the claim that such systems are capable of complex
dynamics and emergent phenomena.

7 Future Work

We are interested in examining the transient growth of recurrent neural networks (RNNs). In the sim-
plest case, we can add a “rounding off ” output layer to discretize the configuration space. Then, IL
transient classification could be used to study the dynamics of RNNs, possibly guiding us toward ap-
propriate network initializations and overall architectures yielding complex dynamics.

It would also be interesting to examine Busy Beaver TMs. These TMs take the longest time to halt
(in the classical sense) among all TMs with the same number of states and tape symbols when run from
an empty tape. It is interesting to observe how such machines behave when run from a randomly
sampled initial configuration and whether they would exhibit complex dynamics, possibly being com-
putationally universal (Zenil, 2012).

Lastly, we could examine the dynamics of systems when simulated from a special region of their
configuration space. The initial configurations could be generated by a designated algorithm, possibly
discovering completely different dynamics of a system, as opposed to its average behavior.

Funding Information
Our work was supported by Grant Schemes at CU, reg. NO. CZ.02.2.69/0.0/0.0/19_073/0016935,
the Ministry of Education, Youth and Sports within the dedicated program ERC CZ under the pro-
ject POSTMAN with reference LL1902, by the Czech project AI&Reasoning CZ.02.1.01/0.0/0.0/
15_003/0000466, European Regional Development Fund, bySVV-2020-260589, and is part of the
RICAIP project that has received funding from the European Unionʼs Horizon 2020 research and
innovation programme under grant agreement No. 857306.

Riferimenti
Aldana, M., Coppersmith, S., & Kadanoff, l. P. (2003). Boolean dynamics with random couplings. In E.

Kaplan, J. E. Marsden, & K. R. Sreenivasan (Eds.), Perspectives and problems in nonlinear science: A celebratory
volume in honor of Lawrence Sirovich (pag. 23–90). Springer. https://doi.org/10.1007/978-0-387-21789-5_2

Bhattacharjya, A., & Liang, S. (1996). Median attractor and transients in random Boolean nets. Physica D:

Nonlinear Phenomena, 95, 29–34. https://doi.org/10.1016/0167-2789(96)00003-6

Channon, UN. (2006). Unbounded evolutionary dynamics in a system of agents that actively process and

transform their environment. Genetic Programming and Evolvable Machines, 7, 253–281. https://doi.org/10.1007
/s10710-006-9009-3

242

Artificial Life Volume 27, Number 3–4

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Cisneros, H., Sivic, J., & Mikolov, T. (2019). Evolving structures in complex systems. Atti del 2019 IEEE
Symposium Series on Computational Intelligence (pag. 230–238). IEEE. https://doi.org/10.1109/SSCI44817.2019
.9002840

Cook, M. (2004). Universality in elementary cellular automata. Complex Systems, 15(1), 1–40.

Crutchfield, J., & Hanson, J. (1993). Turbulent pattern bases for cellular automata. Physica D: Nonlinear

Phenomena, 69, 279–301. https://doi.org/10.1016/0167-2789(93)90092-F

Culik, K., & Yu, S. (1988). Undecidability of CA classification schemes. Complex Systems, 2, 177–190.

Derrida, B., & Flyvbjerg, H. (1987). The random map model: A disordered model with deterministic dynamics.

Journal de Physique, 48, 971–978. https://doi.org/10.1051/jphys:01987004806097100

Derrida, B., & Pomeau, Y. (1986). Random networks of automata: A simple annealed approximation.

Europhysics Letters, 1(2), 45–49. https://doi.org/10.1209/0295-5075/1/2/001

Flyvbjerg, H., & Kjaer, N. J. (1988). Exact solution of Kauffmanʼs model with connectivity one. Journal of Physics

UN: Mathematical and General, 21(7), 1695–1718. https://doi.org/10.1088/0305-4470/21/7/031
Gardener, M. (1970). The fantastic combinations of John Conwayʼs new solitaire game “Life.” Scientific

American, 223, 120–123. https://doi.org/10.1038/scientificamerican1070-120

Gershenson, C. (2004). Introduction to random Boolean networks. In M. Bedau, P. Husbands, T. Hutton, S.
Kumar, & H. Suzuki (Eds.), Proceeding of the workshops and tutorials of the ninth international conference on the
simulation and synthesis of living systems (ALife IX), (Boston, MA, pag. 160–173). CON Premere.

Gutowitz, H. (1990). A hierarchical classification of cellular automata. Physica D: Nonlinear Phenomena, 45(1),

136–156. https://doi.org/10.1016/0167-2789(90)90179-S

Gutowitz, H. (1994). Transients, cycles, and complexity in cellular automata. Physical Review A, 44(12),

R7881–R7884. https://doi.org/10.1103/PhysRevA.44.R7881, PubMed: 9906020

Gutowitz, H. A., Victor, J. D., & Knight, B. W. (1987). Local structure theory for cellular automata. Physica D:

Nonlinear Phenomena, 28(1), 18–48. https://doi.org/10.1016/0167-2789(87)90120-5

Hanson, J. (2009). Cellular automata, Emergent phenomena. In R. UN. Meyers (Ed.), Encyclopedia of Complexity and

Systems Science (pag. 325–335). Springer. https://doi.org/10.1007/978-0-387-30440-3_51

Harris, B. (1960). Probability distributions related to random mappings. Annals of Mathematical Statistics, 31,

1045–1062. https://doi.org/10.1214/aoms/1177705677

Hedlund, G. UN. (1969). Endomorphisms and automorphisms of the shift dynamical system. Mathematical Systems

Theory, 3, 320–375. https://doi.org/10.1007/BF01691062

Huang, S., & Ingber, D. E. (2000). Shape-dependent control of cell growth, differentiation, and apoptosis:
Switching between attractors in cell regulatory networks. Experimental Cell Research, 261(1), 91–103.
https://doi.org/10.1006/excr.2000.5044, PubMed: 11082279

Hudcová, B., & Mikolov, T. (2020). Classification of complex systems based on transients. Negli Atti del
ALIFE 2020: IL 2020 conference on artificial life (Online, Luglio 2020, pag. 367–375). CON Premere. https://doi.org
/10.1162/isal_a_00260

Kaneko, K. (1986). Attractors, basin structures and information processing in cellular automata. In S. Wolfram

(Ed.), Theory and applications of cellular automata (pag. 367–399). World Scientific.

Kari, J. (2005). Theory of cellular automata: A survey. Theoretical Computer Science, 334(1), 3–33. https://doi.org

/10.1016/j.tcs.2004.11.021

Kauffman, S. UN. (1969). Metabolic stability and epigenesis in randomly constructed genetic nets. Journal of

Theoretical Biology, 22(3), 437–467. https://doi.org/10.1016/0022-5193(69)90015-0

Kauffman, S. A., & Weinberger, E. D. (1989). The NK model of rugged fitness landscapes and its application
to maturation of the immune response. Journal of Theoretical Biology, 141(2), 211–245. https://doi.org/10.1016
/S0022-5193(89)80019-0, PubMed: 2632988

Langton, C. G. (1984). Self-reproduction in cellular automata. Physica D: Nonlinear Phenomena, 10(1), 135–144.

https://doi.org/10.1016/0167-2789(84)90256-2

Langton, C. G. (1986). Studying artificial life with cellular automata. Physica D: Nonlinear Phenomena, 22(1–3),

120–149. https://doi.org/10.1016/0167-2789(86)90237-X

Artificial Life Volume 27, Number 3–4

243

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
UN
R
T
l
/

UN
R
T
io
C
e
–
P
D

F
/

2
7
3
–
4
2
2
0
2
0
0
3
3
0
3
UN
R
T
l

_
UN
_
0
0
3
4
2
P
D

B
sì
G
tu
e
S
T

o
N
0
8
S
e
P
e
M
B
e
R
2
0
2
3

B. Hudcová and T. Mikolov

Classification of Discrete Dynamical Systems

Luque, B., & Solé, R. (1997). Phase transitions in random networks: Simple analytic determination of critical

points. Physical Review E, 55, 257–260. https://doi.org/10.1103/PhysRevE.55.257

Luque, B., & Solé, R. (2000). Lyapunov exponents in random Boolean networks. Physica A: Statistical Mechanics

and Its Applications, 284(1), 33–45. https://doi.org/10.1016/S0378-4371(00)00184-9

Lynch, J. F. (1993). A criterion for stability in random Boolean cellular automata. arXiv:adap-org/9305001

Martin, O., Odlyzko, A., & Wolfram, S. (1984). Algebraic properties of cellular automata. Communications in

Mathematical Physics, 93, 219–258. https://doi.org/10.1007/BF01223745

Mitchell, M. (1998). Computation in cellular automata: A selected review. In T. Gramss, S. Bornholdt, M.

Gross, M. Mitchell, & T. Pellizzari (Eds.), Nonstandard computation (pag. 95–140). VCH Verlagsgesellschaft.
https://doi.org/10.1002/3527602968.ch4

Mitchell, M., Crutchfield, J., & Das, R. (2000). Evolving cellular automata with genetic algorithms: A review of
recent work. In Proceedings of the first international conference on evolutionary computation and its applications (EvCAʼ96).
Russian Academy of Sciences.

Neumann, J. V., & Burks, UN. W. (1966). Theory of self-reproducing automata. University of Illinois Press, Urbana, USA.

Ofria, C., & Wilke, C. (2004). Avida: A software platform for research in computational evolutionary biology.

Artificial Life, 10(2), 191–229. https://doi.org/10.1162/106454604773563612, PubMed: 15107231

Owen, UN. B. (2013). Monte Carlo theory, methods and examples. https://statweb.stanford.edu/~owen/mc/

Pineda, O., Kim, H., & Gershenson, C. (2019). A novel antifragility measure based on satisfaction and its

application to random and biological Boolean networks. Complexity, 2019, Article 3728621. https://doi.org
/10.1155/2019/3728621

Ray, T. S. (1991). An approach to the synthesis of life. In C. G. Langton, C. Taylor, J. D. Farmer, & S.

Rasmussen (Eds.), Artificial life II: Proceedings of the workshop on artificial life (Santa Fe, NM, Febbraio 1990,
pag. 371–408). Westview Press.

Reggia, J., Armentrout, S., Chou, H., & Peng, Y. (1993). Simple systems that exhibit self-directed replication.

Scienza, 259(5099), 1282–1287. https://doi.org/10.1126/science.259.5099.1282, PubMed: 17732248

Rogozhin, Y. (1996). Small universal Turing machines. Theoretical Computer Science, 168(2), 215–240. https://doi

.org/10.1016/S0304-3975(96)00077-1

Soare, R. IO. (2016). Turing computability, theory and applications. Springer-Verlag. https://doi.org/10.1007/978-3

-642-31933-4

Soros, l. B., & Stanley, K. O. (2012). Identifying necessary conditions for open-ended evolution through the
artificial life world of Chromaria. In Proceedings of the ALIFE 14: The fourteenth international conference on the
synthesis and simulation of living systems (Manhattan, NY, pag. 793–800). CON Premere.

Stepney, S. (2012). Nonclassical computation: A dynamical systems perspective. In G. Rozenberg, T. Bäck, &
J. N. Kok (Eds.), Handbook of natural computing (pag. 1979–2025). Springer. https://doi.org/10.1007/978-3-540
-92910-9_59

Toffoli, T. (1977). Computation and construction universality of reversible cellular automata. Journal of Computer

and System Sciences, 15(2), 213–231. https://doi.org/10.1016/S0022-0000(77)80007-X

Turing, UN. M. (1937). On computable numbers with as application to the Entscheidungsproblem. Proceedings of

the London Mathematical Society, s2-42(1), 230–265. https://doi.org/10.1112/plms/s2-42.1.230

Vichniac, G. Y. (1984). Simulating physics with cellular automata. Physica D: Nonlinear Phenomena, 10(1), 96–116.

https://doi.org/10.1016/0167-2789(84)90253-7

Wang, X., Lizier, J., & Prokopenko, M. (2011). Fisher information at the edge of chaos in random Boolean
networks. Artificial Life, 17(4), 315–329. https://doi.org/10.1162/artl_a_00041, PubMed: 21762019
Wolfram, S. (1984). Universality and complexity in cellular automata. Physica D: Nonlinear Phenomena, 10(1), 1–35.

https://doi.org/10.1016/0167-2789(84)90245-8

Wolfram, S. (2002). A new kind of science. Wolfram Media.

Wuensche, UN. (1996). The emergence of memory: Categorisation far from equilibrium. In S. R. Hameroff,
UN. W. Kaszniak, & UN. C. Scott (Eds.), Towards a science of consciousness: The first Tucson discussions and debates
(pag. 383–392). CON Premere.