ISSN 13200682 
Source:  http://www.complexity.org.au/ci/vol06/wuensche/wuensche.html  Received:  01/07/1998  
Vol 6:  Copyright 1998  Accepted for publication:  15/10/1998 
Andrew Wuensche
Santa Fe Institute, 1399 Hyde Park Road,
Santa Fe, New Mexico 87501 USA,
Email: wuensch@santafe.edu
WWW: www.santafe.edu/~wuensch/
A key notion in the study of network dynamics is that statespace is connected into basins of attraction . Convergence in attractor basins correlates with order complexity chaos measures on spacetime patterns. A network's ``memory'', its ability to categorize, is provided by the configuration of its separate basins, trees and subtrees. Based on computer simulations using the software Discrete Dynamics Lab [19], this paper provides an overview of recent work describing some of the issues, methods, measures, results, applications and conjectures.
Cellular automata (CA) are a powerful yet simple class of network, characterized by a homogeneous rule and uniform nearest neighbour connections, providing models to study processes in physical systems such as reactiondiffusion [15], and selforganization by the emergence of coherent interacting structures[18]. By contrast, random Boolean networks (RBN) provide models for biological systems such as neural [3] and genetic [11] networks, where connections and rules must be less constrained. In addition, the idealized networks themselves hold intrinsic interest as mathematical/physical systems.
A key notion underlying network behavior is that statespace is organized into a number of basins of attraction, connecting states according to their transitions, and summing up the network's global dynamics, analogous to Poincaré's ``phase portrait'' which provided powerful insights in continuous dynamics.
Figure 1: A basin of attraction (one of 15) of a random Boolean network (n=13, k=3) shown in figure 15. The basin links 604 states, of which 523 are gardenofEden states. The attractor period = 7, and one of the attractor states is shown in detail as a bit pattern. The direction of time is inwards from gardenofEden states to the attractor, then clockwise.
The quality of dynamical behaviour of CA, from ordered to chaotic (1), is reflected by convergence in attractor basins, their characteristic indegree , which influences the length of transients and attractor cycles. The indegree of a state is its number of preimages (predecessors). Bushy subtrees with high indegree imply high convergence and order. Sparsely branching subtrees imply low convergence and chaos. In the case of RBN, attractor basins reveal how the network is able to hierarchically categorize statespace into separate basins, trees and subtrees, the network's ``memory''. Changes to the network's wiring or rules change the memory categories, providing insights into learning[17, 20].
Traditionally, network dynamics has been investigated by running networks forward from many initial states to study spacetime phenomenology[15], and for statistical measures on basins of attraction[8]. More recently, exact representations of basins of attraction and subtrees have become accessible, where algorithms directly compute the preimages of network states, allowing the network to be run ``backwards'' to disclose all possible historical paths[16, 17, 21]. Based on computer simulations using the software Discrete Dynamics Lab (DDLab)[19], this paper provides an overview of network architecture, the characteristics of spacetime patterns, the methods and algorithms for reconstructing basins of attraction, and related parameters, measures, results, applications and conjectures, placing the dynamics along particular trajectories in the context of global dynamics.
1d: k=013. The extra asymmetric cell in even k is on the right. The wiring is shown between two timesteps.
2d: k=213 (k=01 as in 1d). Note that k=6 and k=7 define an effectively triangular grid by changing between odd and even rows. The classical von Neumann and Moore neighbourhoods are indicated.
3d: k=613 (k=05 as in 2d), shown looking up into an axonometric cage.
Figure 2: 1d, 2d and 3d neighbourhood templates defined in DDLab. In 2d and 3d, to maximize symmetry, even k does not include the central target cell.
A CA is a very regular network, sometimes described as an artificial universe with its own physics. Cells take inputs from their nearest (and next nearest) neighbours (local ``wiring'') according to a fixed neighbourhood template, so issues of network geometry and boundary conditions are crucial. The same logical rule is applied everywhere. Figure 2 shows neighbourhood templates for 1d, 2d and 3d as applied in DDLab. An RBN relaxes these constraints, allowing arbitrary ``wiring'' and rules, as in figure 3. The number of input wires available to each cell may also vary. However, an RBN architecture can be biased in countless of ways, described in section 4, to constrain wiring or rules. For example an RBN wiring with a constant rule, or local wiring with mixed rules can be a generalisation of a CA.
The wiring and rules can be tailored to very specific requirements, as in models of neural networks in the cortex [3]. The wiring can be constrained within a fixed distance from each cell, which confers meaning to network geometry and boundary conditions, whereas with completely arbitrary wiring the geometry just provides a convenient way of representing the network. A rule mix can be constrained to sample just a few rules or rules with a particular bias, as in genetic network models[5].
Figure 3: Network wiring. top left: 1d, k=3, the wiring is shown between two timesteps. top right: 2d, k=5. bottom: 3d, k=7. In RBN, cells anywhere in the network are wired back to each position in the ``pseudoneighbourhood''.
Hybrid networks can be constructed by putting an RBN within a CA or viceversa. Networks of networks can be set up, with weak interactions so they perturb each other's dynamics. The functionality for setting up networks in these ways is present in DDLab.
The network parameters can be listed as follows:
A CA neighborhood, or RBN pseudoneighborhood, of size k has permutations of values. The most general expression of the Boolean function or rule is a lookup table (the ruletable) with entries, giving possible rules. Subcategories of rules can also be expressed as simple algorithms, concise AND/OR/NOT logical statements (which could be implemented as combinatorial circuits), totalistic rules[14] or threshold functions.
By convention[14] the rule table is arranged in descending order of the values of neighborhoods, and the resulting bit string converts to a decimal or hexadecimal rule number. For example the k=3 ruletable for rule 30,
The ruletable for other k values are set out in a corresponding way. rules are referred to by their hexadecimal rule numbers. rules are usually referred to by their more familiar decimal rule numbers.
For a given geometry, the behaviour space of CA depends on the size of rulespace, , though rule symmetries effectively reduce this number. For example, the rules in k=3 rulespace reduce to 88 equivalence classes[16]. The behaviour space of RBN is much greater, taking into account possible permutations of wiring and rule schemes, but there are also RBN equivalence classes relating to these permutations[10]. In general, the number of effectively different RBN of size n cannot exceed (see section 9).
Figure 4: Spacetime patterns of a CA (n=24, k=3, rule 90). 24 timesteps from an initial state with a single central 1. Two alternative presentations are shown. Left: cells by value, white=0 black=1. Right: cells colored (or shaded) according to their lookup neighbourhood. This allows filtering, and improves the clarity of spacetime patterns in 2d and 3d.
Figure 5: Spacetime patterns of the k=3 rule 54 (n=150) from the same initial state showing interacting gliders , which are embedded in a complicated background. Left: cells by value. Right: cells by neighbourhood lookup, with the background filtered.
(a) 2d 100x100 triangular grid 
(b) 2d 56x56 +time 
(c) 3d 20x20x20 
Figure 6: Examples of 2d and 3d CA space patterns. (a) is an evolved timestep of a 2d CA on a k=7 triangular lattice with a reactiondiffusion rule. (b) is the 2d gameofLife on a grid, but with a time dimension added, similar to a 1d spacetime pattern. The initial state is set with a number of gliders . (c) is a timestep of a 3d k=7 CA with a randomly selected rule and starting from a single central 1.
Gliders are embedded within a uniform or periodic background or domain, and propagate at various velocities up the system's speed of light set by the neighbourhood diameter. Gliders are interpreted as dislocation in the background or as the boundary reconciling two different backgrounds[4, 18]. Gliders may absorb or eject subgliders (gliderguns). Compound gliders may emerge made up of subgliders recolliding periodically. Figure 3.1 shows some examples.
Glider dynamics has been interpreted as occurring at a phase transition in rulespace between order and chaos [9], relative to the rule parameters [9] and Z[16] (see section 6.2). Inputentropy provides a measure on spacetime dynamics that allows the automatic classification of rulespace (see below).
(a)7e8696de

(b)89ed7106

(c)89ed7106

(d)b51e9ce8

(e)b51e9ce8

Figure 7: Gliders, glider guns and compound gliders in k=5 1d CA . (c) is a compound glider made up of two independent gliders locked into a cycle of repeating collisions. (d) is a glider with a period of 106 timesteps. (e) is a compound glidergun.
Keeping track of the frequency of ruletable lookups (the kblock frequency, or ``lookup frequency'') in a window of timesteps, provides a measure, the variance of inputentropy over time, which is used to classify 1d CA automatically for a spectrum of ordered, complex and chaotic dynamics[22]. The method allows screening out rules that support glider dynamics and related complex rules, giving an unlimited source for further study. The method also shows the distribution of rule classes in the rulespaces of varying neighbourhood sizes. The classification produced seems to correspond to our subjective view of spacetime dynamics, and to global measures on the ``bushiness'' of typical subtrees in attractor basins, characterized by the distribution of indegree sizes in their branching structure.
The lookup frequency can be represented by a histogram (figure 8) which distributes the total of lookups among the neighbourhoods (shown as the fraction of total lookups), where n=system size, w=the window of timesteps defined and k=neighbourhood size. The Shannon entropy of this frequency distribution, the ``inputentropy'' S, at timestep t, for one timestep (w=1), is given by, , where is the lookup frequency of neighbourhood i at time t. In practice the measures are smoothed by being taken over a moving window of timesteps (w=10 in figure 8).
Figure 8: Typical 1d CA Spacetime patterns showing ordered, complex and chaotic dynamics (n=150, k=5, rule numbers shown in hex). Alongside each spacetime pattern is a plot of the inputentropy, where only complex dynamics (centre) exhibits high variance because glider collisions make new gliders.
Figure 8 shows typical examples of ordered, complex and chaotic dynamics in 1d CA, with inputentropy plots and a snapshot of the lookup frequency histogram alongside. In ordered dynamics the entropy quickly settles at a low value with low or zero variance. In chaotic dynamics the entropy settles at a high value, but again with low variance. Both ordered and chaotic dynamics have low inputentropy variance. By contrast, in complex dynamics the entropy fluctuates erratically both up and down for an extended time, because glider collisions produce new gliders, often via a temporary zone of chaotic dynamics. Complex rules can be recognized by eye, subjectively. Inputentropy variance provides a nonsubjective measure for recognizing complex rules automatically.
Figure 9: Entropydensity scatter plot. Inputentropy is plotted against the density of 1s relative to a moving window of timesteps w=10. k=5, n=150. Plots for a number of complex rules from the automatic sample (section 3.3) are show superimposed, each of which has its own distinctive signature, with a marked vertical extent, i.e. high inputentropy variance. About 1000 timesteps are plotted from several random initial states for each rule.
A related method of visualizing the entropyvariance is to plot inputentropy against the density of 1s relative to a moving window of timesteps. Each rule produces a characteristic cloud of points which lie within a parabolic envelope because high entropy is most probable at medium density, low entropy at either low or high density. Each complex rule produces a plot with its own distinctive signature, with high inputentropy variance. Chaotic rules, on the other hand, will give a flat, compact cloud at high entropy (at the top of the parabola). For ordered rules the entropy rapidly falls off with very few data points because the system moves rapidly to an attractor.
Figure 10: top: Classifying a random sample of k=5 rules by plotting mean entropy against standard deviation of the entropy, with the frequency of rules within a grid shown vertically. middle and bottom: Equivalent plots for samples of k=6 and 7 rules.
To distinguish ordered, complex and chaotic rules automatically, the mean inputentropy taken over a span of timesteps is plotted against the standard deviation of the input entropy. Figure 10 summarizes how random samples of k=5, 6 any 7 rules where classified by this method. For each rule, the data was gathered from 5 runs from random initial states, for 430 timesteps, discounting the first 30 to allow the system to settle, with w=5 as the size of the moving window of timesteps.
Chaotic rules are concentrated in the top left corner ``tower", ordered rules on the left with lower entropy. Complex rules have higher standard deviation, and are spread out towards the right. There is a fairly distinct boundary between ordered and chaotic rules, but a gradual transition from both towards the complex rules. As the standard deviation decreases glider interactions either become more frequent, transients longer, tending towards chaos, or less frequent, transients shorter, tending towards order. The plots for k=6 and k=7 rules indicate a greater frequency of chaotic rules at the expense of ordered and complex rules for greater k. The decrease in ordered rules is especially marked.
To check whether the expected dynamics (recognized subjectively) corresponds to the measures as plotted, the dynamics of particular rules at different positions on the plots can be easily examined in DDLab, for example with a mouse click on the scatter plot. Preliminary scans confirm that the expected behaviour is indeed found, but further investigation is required to properly demarcate the space between ordered, complex and chaotic rules and to estimate the proportion of different rule classes for different k.
Input entropy is a local measure on the spacetime patterns of typical trajectories. The distribution of the rule samples according to these local measures may be compared with global measures on convergence in attractor basins, Gdensity and the indegree frequency, described in section 8. Preliminary results indicate a strong relationship between these global measures and the rule sample inputentropy plots.
Figure 11: Spacetime patterns for intermediate 1d architecture, from CA to RBN. n=150, k=5, 150 timesteps from a random initial state. (a) Starting off as a complex CA (rule 6c1e53a8 as in figure 8), 4% (30/750) of available wires are randomized at 30 timestep intervals. The coherent pattern is progressively degraded. (b) A network with local wiring but mixed rules, vertical features are evident. (c) RBN, random wiring and mixed rules, with no bias, shows maximal chaotic dynamics.
In contrast to CA, glider dynamics in general cannot occur in RBN because of their irregular architecture. Figure 11 (left) shows glider dynamics degrading as local wiring is progressively scrambled. An alternative orderchaos notion in RBN is the balance between ``frozen'', stabilized, regions and changing regions in the spacetime pattern[8]. Stable regions are characteristic of RBN with low connectivity, , because rules which induce stability are relatively frequent in these rulespaces. To induce stability for , where chaotic rules become overwhelmingly predominant, biases on rules must be imposed, low (see section 6.2) or a high proportion of ``canalizing'' inputs. In a rule's lookup table, an input wire is canalizing if a particular input (0 or 1) determines, by itself, the neighbourhood's output. A rule's degree of canalization can be from 0 to k, for the same output; for the network it is the percentage of all inputs that are canalizing, C. An RBN's orderchaos characteristics, for varying C, are captured by the measures illustrated in figure 12, and described below.
Figure 12: Orderchaos measures for a RBN , k=5. C = the percentage of canalizing inputs in the randomly biased network. top left: frozen elements that have stabilized for 20 timesteps are shown, 0sgreen, 1s red, otherwise white, for C=25% and 52%. top right: the loglog ``damage spread'' histogram for C=52%, sample size about 1000. bottom: the Derrida plot for C=0%, 25%, 52%, and 75%, for 1 timestep, =00.3, interval = 5, sample for each = 25.
The ``Derrida plot''[2] , is analogous to the Lyapunov exponent in continuous dynamics, and measures the divergence of trajectories based on normalized Hamming distance H, the fraction of bits that differ between two patterns. Pairs of random states separated by , are independently iterated forward by one (or more) timesteps. For a sample of random pairs, the average is plotted against , and the plot is repeated for increasing (from 0 to 0.3 in figure 12). A curve above the main diagonal indicates divergent trajectories and chaos, below  convergence and order. A curve tangential to the main diagonal indicates a balance. A related measure is the distribution of ``damage spread'' resulting from a single bit change at a random position in a random state, for a sample of random states. The size of damage is measured once it has stabilized, i.e. not changed for say 5 timesteps. A histogram (figure 12) is plotted of damage size against the frequency of sizes. Its shape indicates order or chaos in the network, where a balance between order and chaos approximates to a power law distribution. Results by these measures for k=5, indicate a balance at (see figure 12). There are further measures on basins of attraction as in figure 14.
These methods are applied in the context of RBN models of genetic regulatory networks[8] discussed in section 10. The conjecture is that evolution maintains genetic regulatory networks marginally on the ordered side of the orderchaos boundary to achieve stability and adaptability in the pattern of gene expression which defines the cell type[5].
For a network size n, an example of one of its states B might be . Statespace is made up of all states, the space of all possible bitstrings or patterns.  
Part of a trajectory in statespace, where C is a successor of B, and A is a preimage of B, according to the dynamics of the network.  
The state B may have other preimages besides A, the total is the indegree. The preimage states may have their own preimages or none. States without preimages are known as gardenofEden states.  
Any trajectory must sooner or later encounter a state that occurred previously  it has entered an attractor cycle. The trajectory leading to the attractor a transient. The period of the attractor is the number of states in its cycle, which may be only just one  a point attractor.  
Take a state on the attractor, find its preimages (excluding the preimage on the attractor). Now find the preimages of each preimage, and so on, until all gardenofEden states are reached. The graph of linked states is a transient tree rooted on the attractor state. Part of the transient tree is a subtree defined by its root.  
Construct each transient tree (if any) from each attractor state. The complete graph is the basin of attraction. Some basins of attraction have no transient trees, just the bare ``attractor''.  
Now find every attractor cycle in statespace and construct its basin of attraction. This is the basin of attraction field containing all states in statespace, but now linked according to the dynamics of the network. Each discrete dynamical network imposes a particular basin of attraction field on statespace. 
Figure 13: State space and basins of attraction.
The idea of basins of attraction in discrete dynamical networks is summarized in figure 13. Given invariant network architecture and the absence of noise, a discrete dynamical network is deterministic, and follows a unique (though in general, unpredictable) trajectory from any initial state. When a state that occurred previously is revisited, which must happen in a finite statespace, the dynamics becomes trapped in a perpetual cycle of repetitions defining the attractor (state cycle) and its period (minimum one, a stable point).
These systems are dissipative. A state may have multiple ``preimages'' (predecessors), or none, but just one successor. The number of preimages is the state's ``indegree''. Indegrees greater than one require that transient states exist outside the attractor. Tracing connections backwards to successive preimages of transient states will reveals a treelike topology where the ``leaves'' are states without preimages, known as gardenofEden states. Conversely, the flow in statespace is convergent. Measures of convergence are Gdensity, the fraction of states that are gardenofEden, and the distribution of indegrees, described in section 8. The set of transient trees rooted on the attractor is its basin of attraction (figure 1). The local dynamics connects statespace into a number of basins, the basin of attraction field, representing the systems global dynamics (figure 15).
Figure 14: Statistical data on attractor basins for a large network; a 2d RBN 20 20, k=5, with fully random wiring and a fraction of canalizing inputs C=50%. The histogram shows attractor types and the frequency of reaching each type from 12,360 random initial states, sorted by frequency. 46 different attractors types where found, their periods ranging from 4 to 102, with average transient length from 21 to 113 timesteps. The frequency of arriving at each attractor type indicates the relative size of the basin of attraction.
Attractor basins are constructed with algorithms that directly compute the preimages of network states[16, 17, 21]. This allows the network's dynamics, in effect, to be run backwards in time. Backward trajectories will, as a rule, diverge. Different reverse algorithms apply to networks with different sorts of connectivity. The most computationally efficient algorithm applies to 1d networks with local wiring, taking advantage of the regularity of connections. The wiring must be uniform, as for 1d CA, but the network may have a mix of rules. Analogous algorithms could be derived for 2d and 3d networks, but have not been implemented. An alternative algorithm is required for RBN with their nonlocal connections and possibly mixed k. This algorithm also applies to CA of any dimension or geometry, as CA are just a subclass of RBN.
Provided , these methods are in general orders of magnitude faster than the brute force method (section 9.1), constructing an exhaustive map resulting from network dynamics, a method which rapidly becomes intractable with increasing network size and so is limited to very small systems. However, the exhaustive method may be applied to all types of network, and also allows the attractor basins of random maps to be constructed, as described in section 9. The agreement of these three independent methods, and other checks, give considerable confidence in the accuracy of the preimage computations.
Some basic information on attractor basin structure can be found by statistical methods, first applied by Walker[13], as shown in figure 6. These are also implemented in DDLab and are appropriate for large networks. Trajectories are run forward from many random initial states looking for a repeat in the network pattern to identify the range of attractor types reached. The frequency of reaching a given attractor type indicates the relative size of the basin of attraction, and other data are extracted such as the number of basins, and the length of transients and attractor cycles.
Figure 15: The basin of attraction field of a random Boolean network (n=13, k=3). The states in state space are organized into 15 basins, with attractor periods ranging between 1 and 7. The number of states in each basin is: 68, 984, 784, 1300, 264, 76, 316, 120, 64, 120, 256, 2724, 604, 84, 428. Figure 1 shows the arrowed basin in more detail. Right: the network's architecture, its wiring/rule scheme. 
Consider a 1d CA of size n (indexed ) and neighbourhood k. To find the all preimages of a state A, let P be a ``partial preimage'' where at least k1 continuous bits (on the left) up to and including , are known. Now find the next unknown bit to the right, , consistent with the ruletable. ( indicates known, unknown, bits),
If k=3 (for example), the bitstring corresponds to two neighbourhood entries in the ruletable. When their outputs, and , are compared with each other and with there are three possible consequences. The permutation is either deterministic, ambiguous or forbidden.
If forbidden (3) the partial preimage P is rejected. If deterministic or ambiguous (1 or 2) the procedure is continued to find the next unknown bit to the right. However, in the ambiguous case (2), both alternative partial preimages must be continued. In practice one is assigned to a stack of partial preimages to be continued at a later stage. As the procedure is reapplied to determine each successive unknown bit towards the right, each incidence of ambiguous permutations will require another partial preimage to be added to the stack. Various refinements can limit this growth.
The procedure is continued to the right to overlap the assumed start string, to check if periodic boundary conditions are satisfied; if so the the preimage is valid. The procedure is reapplied to each partial preimage taken from the partial preimage stack, starting at the first unknown cell. Each time an ambiguous permutation (2) occurs a new partial preimage must be added to the stack, but the stack will eventually be exhausted, at which point all the valid preimages containing the assumed start string will have been found. The procedure is applied for start strings, assuming the different possible values of the first k1 bits. The reverse algorithm is applied from left to right in DDLab, but is equally valid when applied from right to left. Examples are given in [16, 21].
A by product of the CA reverse algorithm is the probability of the next unknown bit being deterministic (section 6.1(1)). Two versions of this probability are calculated from the ruletable. for the reverse algorithm applied from left to right, and for the converse. The Z parameter is the greater of these values. For Z=1 it can be shown[16] that for any system size n, the maximum indegree, , because the next unknown bit is always uniquely determined, so the assumed start string of length k1 may generate at most preimages. If only one of or =1, , because at least one assumed start string must be forbidden (section 6.1(3)). At the other extreme, for Z=0, all state space converges on the state all0s or all1s in one step. For high Z, low indegree (relative to system size n) is expected in attractor basins, growing at a slow rate with respect to n. Conversely, for low Z, high relative indegree is expected growing quickly with respect to n. High Z predicts low convergence and chaos, low Z predicts high convergence and order.
The neighborhoods of size k, each indexed , each have an output T (0 or 1) which makes up the ruletable (section 2), and may be expressed as . To calculate from the rule table, let be the count of ruletable entries belonging to deterministic pairs, such that,
and (not T)
The probability that the next bit is determined because of the above is given by, . This is a first approximation of .
Let be the count of ruletable entries belonging to deterministic 4tuples (where `` '' may be 0 or 1), such that,
and
The probability that the next bit is determined because of the above is given by, . This count is repeated if necessary for deterministic 8tuples where , 16tuples, 32tuples, up to the special case of just one tuple which occupies the whole ruletable. These are are independent nonexclusive probabilities that the next bit is determined. The union of the probabilities , is given by the following expression (the order of the probabilities makes no difference to the result),
which simplifies to,
and may be expressed as(2) where , and = the count of ruletable entries belonging to deterministic tuples. A converse procedure gives , and the Z parameter = the greater of and . Examples are given in [16, 21].
By virtue of being a convergence parameter, Z is also an orderchaos parameter varying from 0(order)  1(chaos). Z can be compared with Langton's[9] well known parameter (3). is an orderchaos parameter for CA which may have values greater than binary, and measures the density of ``nonquiescent'' outputs in a ruletable, so for just binary CA, where c=the count of 1s a ruletable on k inputs. varies between 0 (order)  0.5 (chaos)  1 (order). To allow Z and to be compared, a normalized version of binary is defined[16], where is the count of 0s or 1s in the ruletable, whichever is less. then varies from 0 (order)  1 (chaos) just as Z.
Plots of Gdensity against both the and Z parameters, showing the discrepancies as well as similarities, are shown in figure 6.2, for the 256 k=7 totalistic rules, which reduce to 136 nonequivalent rules in 72 clusters (having equal and Z). Points plotted in the top right corner of the graph represent values that do not correspond to behaviour as expected.
Figure 16: Gdensity against both and Z for the set of k=7 totalistic rules, n=16, for . The complete basin of attraction field was generated for each rule and gardenofEden states counted.
Figure 17: Computing RBN preimages. The changing size of a typical partial preimage stack at successive elements. n=24, k=3.
Consider an RBN of size n. Find all preimages of a state A, ). Each network element , has a pseudoneighbourhood size (assuming a mixed k network), indexed , a wiring scheme , ), where is a number between n1 and 0, the position of the wire connection from the jth branch of the pseudoneighbourhood, and a ruletable .
To find the all preimages of A, let P, , be a candidate preimage consisting of empty network elements, as yet unassigned to either 0 or 1. Starting with an element of A, , assign bits from all valid pseudoneighbourhoods in the ruletable , i.e. that are consistent with , to separate copies of P according to the wiring scheme . As there will be a mix of 0s and 1s in , only some of the possible pseudoneighbourhoods will be valid. This will produce a stack of ``partial preimages'' with some bits allocated and the remainder empty.
Now repeat the procedure for another element of A, say , but this time independently for each partial preimage previously created. If the allocation of a bit to a given partial preimage conflicts with the bit already assigned, then the partial preimage is rejected. Otherwise, the partial preimage is added to the next generation of partial preimages in a new stack. The allocation will be valid if it is made to an ``empty'' element, or to an allocated element with an equal bit. Valid allocation increases the size of the partial preimage stack, conflicts reduce the size of the stack.
This procedure is repeated in turn for the remaining network elements of A. If the stack size is reduced to zero at any stage A has no preimages. The algorithm works for any ordering of elements in A, though to minimizes the growth of the partial preimage stack, the order should correspond to the greatest overlap of wiring schemes. The changing size of the stack at successive elements can be displayed in DDLab, an example is shown in figure 17. When the procedure is complete, the final preimage stack may still have empty network elements, which did not figure in any wiring scheme. These are duplicated so that all possible configurations at empty element positions are represented. The resulting preimage stack is the complete set of preimages of A without duplication.
The reverse algorithm for RBN works for networks with any degree of intermediate architecture between RBN and CA, including CA of any dimension. More detailed explanations of the algorithm are given in [17, 21].
To construct a basin of attraction containing a particular state, the network is iterated forward from the state until a repeat is found and the attractor identified. The transient tree (if it exists) rooted on each attractor state is constructed in turn. Using one of the reverse algorithms, the preimages of the attractor state are computed, ignoring the preimage lying on the attractor itself. Then the preimages of preimages are computed, until all ``gardenofEden'' states have been reached.
In a similar way, just a subtree may be constructed rooted on a state. Because a state chosen at random is very likely to be a gardenofEden state, it is usually necessary to run the network forward by at least one timestep, and use the state reached as the subtree root. Running forward by more steps will reach a state deeper in the subtree so allow a larger subtree to be constructed.
For CA, a considerable speedup in computation is achieved by taking advantage of ``rotational symmetry''[16], a property of the regularity of CA and periodic boundary conditions, resulting in equivalent subtrees and basins.
Attractor basins are portrayed as state transition graphs, vertices (nodes) connected by directed edges. States are represented by nodes, by a bit pattern in 1d or 2d, or as the decimal or hex value of the state. In the graphic convention[16, 19], the length of edges decreases with distance away from the attractor, and the diameter of the attractor cycle approaches an upper limit with increasing period. The direction of edges (i.e. time) is inward from gardenofEden states to the attractor, and then clockwise around the attractor cycle, as shown in figure 1. Typically, the vast majority of states in a basin of attraction lie on transient trees outside the attractor, and the vast majority of these states are gardenofEden states.
Figure 18: The Gdensity plotted against system size n, for the ordered, complex and chaotic rules shown in figures 8 and 19. The the entire basin of attraction field was plotted for n=7 to 22, and gardenofEden states counted. The relative Gdensity and rate of increase with n provides a simple measure of convergence.
Measures on attractor basins include the number of attractors, attractor periods, size of basins, characteristic length of transients and the characteristic branching within trees. The last in particular gives a good measure of the convergence of the dynamical flow in statespace, where high convergence indicates ordered, and low convergence indicates chaotic dynamics.
The simplest measure that captures the degree of convergence is the density of gardenofEden states[18], Gdensity, counted in attractor basins or subtrees, and the rate of increase of Gdensity with n as shown in figure 18. A more comprehensive measure is the indegree frequency distribution, plotted as a histogram. The indegree of a state is the number of its immediate preimages. This can be taken on a basin of attraction field, a single basin, a subtrees, or on just part of a subtree for larger systems. Subtrees are portrayed as graphs showing trajectories merging onto the subtree root state.
Examples of indegree histograms for typical subtrees of ordered, complex, and chaotic rules are shown in figure 19. The horizontal axis represents indegree size, from zero (gardenofEden states) upwards, the vertical axis represents the frequency of the different indegrees. The system size n=50 for the complex and chaotic rules. For very ordered rules indegrees become astronomical. The ordered rule shown is only moderately ordered, however the system size was reduced to n=40 to allow easier computation.
From the preliminary data gathered so far, the profile of the indegree histogram for different classes of rule is as follows:
Ordered rules: Very high gardenofEden frequency and significant frequency of high indegrees. High convergence.
Complex rules: Approximates a power law distribution. Medium convergence.
Chaotic rules: Lower gardenofEden frequency compared to complex rules, and a higher frequency of low in degrees. Low convergence.
Figure 19a: Ordered dynamics. Rule 01dc3610, n=40, Z=0.5625, =0.668. right: The complete subtree 7 levels deep, with 58153 nodes, Gdensity=0.931.  
Figure 19b: Complex dynamics. Rule 6c1e53a8, n=50, Z=0.727, =0.938. right: The subtree, stopped after 12 levels, with 144876 nodes, Gdensity=0.692.  
Figure 19c: Chaotic dynamics. Rule 994a6a65, n=50, Z=0.938, =0.938. right: The subtree, stopped after about 75 levels, with 9446 nodes, Gdensity=0.487. 
The attractor basins of discrete dynamical networks can be put into the wider context of random graph theory. CA belong to the set of RBN which in turn belong to the set of random directed graphs with out degree one, known as random maps. This is a mapping of the Boolean hypercube of sequences of length n (i.e. a set of size comprising all binary strings of length n), a mapping from . The structures found in random maps correspond to those in the attractor basins of discrete dynamical networks, where each separate component of the graph is made up of trees rooted on just one closed cycle. The structures can be computed in DDLab just as the attractor basins of CA or RBN.
Figure 20:The basin of attraction field of a typical unbiased random map, n=12. The states in state space are connected into 9 basins of attraction. The period (p) and size (s) of the biggest three (top row), including their percentage of statespace, are as follows: (1) p=118 s=3204=78.2%. (2) p=20 s=599=14.6%. (3) p=32 s=194=4.74%. The field's Gdensity=0.37, this is a low value implying chaotic dynamics.
A random map can be constructed by assigning a successor to each state in statespace, i.e. independently assign one successor (or image) also belonging to , chosen at random (or with some bias) to each element of the set . There are possible mappings. The mapping is represented below as pairs of strings (or states in statespace), where each image represents a possibly different member of the set .
The list of images is likely to contain repeats, and if so some other members of must be missing from the list. Transitions to some arbitrary element may thus be onetoone or manytoone, or may not exist. The latter is a gardenofEden state in the terminology of discrete dynamical networks. A representation of the particular mapping may be drawn as a basin of attraction field or fragment thereof just as for CA or RBN, and will have the same general topology of trees rooted on attractor cycles, as shown in figure 20.
Random maps provide the most general context for a discrete dynamical system and are equivalent to a fully connected RBN where k=n, (the neighbourhood = the network size). This follows because each cell in the RBN can be assigned an arbitrary output for any network state.
The ``brute force'' reverse algorithm for finding the preimages of states in random maps can also be applied to discrete dynamical networks, RBN and CA. The method depends on first constructing an exhaustive mapping . For discrete dynamical networks, the mapping is defined by iterating the network forward by one step from every state in statespace and filling in the image list accordingly. A list of pairs, each state and its image (successor), is held in a data structure. The preimages of an arbitrary state S are found by scanning the image list; any occurrence of S in the list gives a preimage, the state paired with S. If S does not occur in the list it has no preimages, a gardenofEden state.
Genetic regulatory networks have been thought of as discrete dynamical networks, to explain how gene expression is able to settle into a number of distinct stable patterns or cell types, despite the fact that all eukaryotic cells in an organism carry an identical set of genes[5, 7, 11, 23]. The gene expression pattern of a cell needs to be stable but also adaptable. Section 4 described biases to RBN to achieve such a balance, and related measures.
Cell types have been interpreted as the separate attractors or basins of attraction into which network dynamics settles from various initial states. Trajectories leading to attractors are seen as the pathways of differentiation. The attractor basins in RBN are idealized models for the stability of cell types against mutations, and also perturbations of the current state of gene activation. Figure 21 illustrates both effects. If a particular reference state (pattern of gene activation) undergoes a 1 bit perturbation, the dynamics may return to the same subtree, the same basin, or it may be flipped to another basin, a different cell type. In this case the basin of attraction field remains unchanged. Alternatively, the network itself my undergo a mutation (in the genotype), resulting in an an altered basin of attraction field (the phenotype).
The examples in figure 21 are small so that the pattern at each node can be shown. Larger networks are affected in analogous ways. The consequences of a one bit mutation has a relatively smaller effect with increasing network size. However, a particular one bit mutation may cause drastic consequences whatever the size, such as breaking an attractor cycle. The consequences of moving a connection wire is usually greater than a one bit mutation in a rule.
Attractors classify statespace into broad categories, the network's ``content addressable'' memory in the sense of Hopfield[6]. Furthermore, statespace is categorized along transients, by the root of each subtree forming a hierarchy of subcategories. This notion of memory far from the equilibrium condition of attractors greatly extends the classical concept of memory by attractors alone[17, 20].
It can be argued that in biological networks such as neural networks in the brain or networks of genes regulating the differentiation and adaptive behavior of cells, attractor basins and subtrees, the network's memory, must be just right for effective categorization. The dynamics need to be sufficiently versatile for adaptive behavior but short of chaotic to ensure reliable behavior, and this in turn implies a balance between order and chaos in the network.
A current research topic, known as the ``inverse problem'', is to find ways to deduce network architecture from usually incomplete data on transitions, such as a trajectory. This is significant in genetics, to infer the genetic regulatory network (modeled as RBN) from data on successive patterns of gene expression in the developing embryo[12]. In pattern recognition and similar applications in the area of artificial neural networks, solutions to the inverse problem would provide ``learning'' methods for RBN to make useful categories[17, 20].
Figure 21:The basin of attraction field of (a) The RBN (n=6, k=3) as defined in the table (above), and (b) the RBN following a 1 bit mutation to one of its rules. Some differences in the fields are evident. The result or a 1 bit perturbation to a reference state of all 1s (rs) is indicated by its 1 bit mutants (m). 
Important insights may be gained by considering network dynamics in the
context of attractor basins. Some methods of achieving this have been
presented, including parameters and measures on particular trajectories
that may be related to those on global dynamics. It is hoped that these
methods may provide a basis for future research, both in theory and
applications, in the many areas of complex systems where network
dynamics plays a central role.
Acknowledgments: Thanks to Cosma Shalizi for suggestions, and to the many people who have contributed to this work over the years, notably Chris Langton and Stuart Kauffman.