CONTROLLED CELLULAR AUTOMATA

. Cellular Automata have been successfully used to model evolution of complex systems based on simples rules. In this paper we introduce controlled cellular automata to depict the dynamics of systems with controls that can aﬀect their evolution. Using theory from discrete control systems, we derive results for the control of cellular automata in speciﬁc cases. The paper is mostly oriented toward two applications: ﬁre spreading; morphogenesis and tumor growth. In both cases, we illustrate the impact of a control on the evolution of the system. For the ﬁre, the control is assumed to be either ﬁrelines or ﬁrebreaks to prevent spreading or dumping of water, ﬁre retardant and chemi-cals (foam) on the ﬁre to neutralize it. In the case of cellular growth, the control describes mechanisms used to regulate growth factors and morphogenic events based on the existence of extracellular matrix structures called fractones. The hypothesis is that fractone distribution may coordinate the timing and location of neural cell proliferation, thereby guiding morphogenesis, at several stages of early brain development.


1.
Introduction. The vast range of applications of cellular automata spans numerous disciplines -mathematics, computer science, computer technology, biology, business and many more. Several individual cellular automata are even known to be Turing complete, meaning that any algorithm can be implemented by specifying the initial configuration. This fact alone demonstrates that cellular automata are capable of generating systems of boundless complexity.
Rather than analyzing any single application of cellular automata in detail, or even focusing on a specific set of rules, we consider how we can apply control to cellular automata to guide them toward particular outcomes. We consider controls that fit into two categories: modifications to the states of the cells that do not conform to the rules of the automaton and modifications to the rules that define the automaton being run. Techniques from discrete control systems are used for preliminary results and are illustrated by various examples and simulations. This paper is a first approach to the development of systematic tools to analyze controlled cellular automata. In particular, in Section 3.2, we establish fundamental controlabllity results for cellular automata whose rules are linear which serve as a foundation for future work in the field.
In order to illustrate the potential applications for control, we consider a simple model for the spreading of fire and a model for cell proliferation. Our main motivation comes from the fact that hybrid automata and cellular automata have been increasingly used to mathematically model biological system, indeed CAs are especially well-suited to create complex evolution dynamics with simple rules. For instance, there is significant literature of cellular automata used in tumor growth and related therapy ( [1] and references therein). The model for cell proliferation presented here represents a different approach to the co-evolving cellular automata that the first and second authors developed in [2]. In that conference paper, we introduced a morphogenetic cellular automaton to study how an embryo organizes and directs cellular growth. It is based on the recent discovery of specialized extracellular matrix structures. More precisely, some factors influencing growth have been characterized, in particular forces between cells [4,13] and growth factors detected by the cells through receptors on the cell surface [14]. However, specialized extra-cellular matrix structures, known as fractones [10], have been identified as a control mechanism orchestrating the coordinated diffusion of the growth factors through the extra-cellular space [5,8,11]. They were discovered in the neurogenic zone of the adult mammalian brain [10,12]. Fractones can promote cell proliferation and can also inhibit cell proliferation in the adult brain. More precisely, since the basal laminae is thought to contribute to morphogenesis, it was hypothesized that fractones capture and store growth factors until this concentration reaches a prescribed threshold to initiate binding with cell-surface receptors. It is believed that several types of fractones may exist, each targeting specific growth factors, which are bound, stored, concentrated, and presented to nearby cells, thus triggering a reaction according to the type of growth factor. The distribution of the fractones is model as the control function in the set of rules, and we also add cell death in the model to allow for more complex and refined shapes. The model is then used as a base for tumor growth where two types of cells are considered (STEM cells and tumor cells). Section 4.2 describes the cell and tumor growth applications in details and provides numerical simulations.
2. Cellular automaton. A Cellular Automaton (CA) is a collection of computing cells which repeatedly update their internal states. A CA is defined on a grid G that we will call the cellular space. In this paper, we consider continuous automaton. More precisely, we have the following definition see [2,9,7,6] for discrete version of it. 1. S is the set of possible states for a cell. We consider a state to be string of continuous values, and we assume k to be the length of the string. We have S ⊆ R k and this subset is determined by the state constraints. 2. The constant d is the dimension of the automaton, d = 1, 2 or 3. 3. N defines the neighborhood within which the local interaction of agents occurs. N : G → P(G). The most common ones are the Von Neumann neighborhood and the Moore neighborhood. Given a cell c ∈ G, we denote by n(c) the dimension of its neighborhood (typically n(c) is constant over G but not necessarily). 4. f is a function specifying the transition rule governing changes in an agent's state. The output depends on the state of the agent and its neighbors.
Note that in our definition the cells remain discretely separated from each other and we will consider a discrete time evolution. Most common in the literature is a definition with S = [0, 1] where the state value is a bounded single value, but we are extending to multi-valued states and possibly with no state constraints.
We introduce s : G → S G the state function which associates to a cell c its states, and we denote by s t (c) the state of cell c at time t. The function s defines the notion of configuration for our automaton.
2.1. Metric. To compare two configurations for a given CA, we need a definition of distance between them. They are multiple ways to do that, and the distance we present is based on our intuition of what we think is relevant to our applications.
2. An extended Moore neighborhood of radius r, denoted N r M , is a 2r + 1 × 2r + 1 square-shaped neighborhood which contains (2r + 1) 2 + 1 cells. Note that in our definition, the cell c is contained in the neighborhood. Definition 2.3. Let s 1 t and s 2 t be two configurations at time t of a given CA. We introduce with p = 1, · · · , k be the component p of the state's string. The distance between s 1 t and s 2 t is given by: where P (c) is the cardinality of the set of neighbors of c over the three neighborhood N r M (c). It will be typically 35 = 2 r=0 (2r + 1) 2 unless it is close to the boundary of the grid G and will therefore be smaller. Proof. By definition d is non-negative, and is symmetric: Example 1. To illustrate our metric, we use a simple example associated to a m × m grid with single binary state s(c) ∈ {0, 1}. The configuration s 0 to which we will measure the distances from is the classical checkerboard, see Fig. 1 left. If we denote by c ij the cell in unit i, j of the grid with i representing the rows and j the columns, we have that: We associate the color black to the number 1 and the color white to 0. We now consider the checkerboard such that the colors are alternated compare to s 0 , see Fig. 1 right, i.e.: Applying the formula, we can write explicitly the formula which can be approximated by: where we assume for simplicity of the formula the cardinality of the set of neighbors to be P (c) = 35 for all cells even the ones on the boundary. Let us consider an additional configuration s 2 with all cell having the same value, for instance we assume s 2 (c ij ) = 0 for all i, j ∈ {1, · · · , m} which means that all cells are white. It can be shown that d(s 0 , s 2 ) ≈ 1 2 and therefore d(s 0 , s 1 ) ≈ 3 35 = 0.086 is much smaller than d(s 0 , s 2 ) = 0.5. The interpretation is as follows.
Even though configurations s 0 and s 1 are such that their value differs in each cell, their appearance is the same (both are checkerboard with alternating colors). We therefore want to consider to be close configurations and our choice of distance reflects this. On the other side configurations s 0 and s 2 differ only for half of the cells, however their appearance is much different since one is a checkerboard and the other is uni-color. Our distance does take this into account and measure s 0 and s 2 as much further apart than s 0 and s 1 . This shows us that our metric depends on the structure of automata more than it depends on the particular values in cells.
Definition 2.5. Let s 1 and s 2 be two arbitrary configurations of a given CA. We introduce the difference between them as (s 1 − s 2 ) p (c) = s 1 p (c) − s 2 p (c). Note that s 1 − s 2 might not be a configuration of CA since it does not necessarily satisfy the constraint that s 1 (c) − s 2 (c) ∈ S.
Note that the distance defined in 2.3 can be applied to any two configurations even if they do not belong to the same CA. It therefore make sense to compute d(s 1 − s 2 ), s 1 , s 2 ∈ CA even though s 1 − s 2 might not belong to CA. Clearly, we have d(s 1 , s 2 ) = d(s 1 − s 2 , 0) where 0 is the configuration taking the value 0 in every cell.
The next lemma illustrates that if the difference in state's values of two configurations are clustered, these two configurations are considered closer together than if the same difference in state's values was more spread out over the grid. See illustration in Figure 2. We introduce the equivalence relation ∼ as follows. We say that s 1 − s 0 ∼ s 2 − s 0 if for each value associated to a cell in s 1 − s 0 , there exists a

Controlled transition function.
Definition 2.7. The uncontrolled transition function is defined as: where n associates to each cell c the dimension of its neighborhood. The global dynamics is given by: By definition the transition function is a map whose argument is the states of the cells in the neighborhood of a prescribed central cell c and image is the new state of the cell c: The function f can take many forms depending on the application under study.
Typical states constraints are of the form c p,min ≤ s p (c) ≤ c p,max , c p,min ≤ 0 ≤ c p,max for all c ∈ G and each component of the state string p = 1, · · · , k. In this case, the transition between two states is expressed as: where α p i,c : I → R are real valued functions defined on I the discrete time interval. Again, in case of state constraints the transition between states need to be modified accordingly.
In our work we consider an external parameter that can affect the transition function and therefore the evolution of the system. Definition 2.8. The control of a cellular automaton is defined as a discrete-time map: where l depends on the control mechanism for the application under study, see Section 4 for specific examples. The control domain accounts for practical constraints on the control, such as a bound on its magnitude. With presence of a control, the transition function becomes: and the global dynamic is now F (s t , u t )(c) = f (s t (N (c), u t (c))) = s t+1 (c).

Evolution and control problem.
3.1. Evolution. We have to define the evolution (trajectory) for our system. Let I = {0, 1, · · · , T } be a discrete time interval, and an initial configuration s 0 : G → S G . The evolution Evol s0,u T for a given control u is a sequence of configuration s t , t = 0 · · · , T given by the global dynamics: where F t represents t composition of the map F .
Note that for a cell on the border of the grid the number of neighbors will be smaller. Assume that the coefficients do not vary with time and that p = 1, we can therefore eliminate the dependence on p in the transit function which is now expressed as: Assume the following coefficients: α 00 = 0, α 10 = 0.24, α −10 = 0.24, α 01 = 0.24 and α 0−1 = 0.24. We chose the initial configuration with ones in the diagonal and zeroes elsewhere: Figure 3 illustrates the evolution of this CA at time t = 0, t = 5, t = 20 and t = 100. Notice the spreading and exponential decay that occurs during the simulation. Assume we have a control that can alter the transition rule over the entire grid. The control is designed to be an on-off switch, u(t) ∈ {0, 1} and is implemented as follows: The control changes the weight at which a cell's evolution is impacted by the cells in its neigborhood. Since it is a scalar control there are only two possible transition rules though additional controls can be added to offer more flexibility as well as the possibility to impact only a subset of the cells in the grid. In the sequel we will therefore focus on situations where the control acts in specific ways. This is the case when for instance: The control switches on at 20, off at 40, back on at 50 and finally off again at 100. Notice that while the control is off, the sign is constant and the magnitude diminishes at a slow exponential rate. When the control is on, the sign alternates with each timestep and the magnitude increases at a slow exponential rate.
• the control can act only on a subset of the set of cells in the grid; • the control is bounded; • the value a control can take on one component of a state is constant over all cells in the grid; • etc.
For instance, the control can act on the grid as: where B t is an k × l matrix and u t is a l × 1 vector. The matrix B t indicates how the control is acting on the k components of the state of cell c (B t ≡ 0 means that no control is applied to c at time t) and u t takes its value in the control domain (determined by the application) which is a subset of R l .  To find a general framework to study the controllability problem, we will express our cellular automaton as a discrete control system. We introduce a configuration vector q that represents the value of each component of the state of a cell for all cells in the square lattice G. Since the state of a cell has k components and we assume the grid is m × m, our state vector q belongs to a subset (R m × R m ) k which is determined by the state constraints. We introduce n = m 2 k the dimension of vector q. The global dynamics can then be written as: 3.2.1. Linear case. Assume first that the control is fixed to 0 (u ≡ 0), the evolution of the system is governed by the uncontrolled transition function (13). Our system can be written as a linear discrete-time system: where A t is a n × n matrix, B t is a n × l matrix and l the dimension of the control (typically a multiple of k as in section 4.1). The A matrix is typically a sparse matrix since the non zero elements occur only through the notion of neighbors, and in most problems the rules of the uncontrolled transition function do not depend on time, which makes A constant. When both A and B are constants, the system is said to be time invariant. Clearly, we have the following result.
Proposition 1. The solution to (23) is given by: where the state transition matrix φ(t, j) is defined by Time Invariant System. In this case the solution is given by where A i denotes the matrix A to the power i. Assume that A is diagonalizable, which means that there is an invertible matrix P such that A = P DP −1 and D is diagonal. A necessary and sufficient condition is that A has n linearly independent eigenvectors which is true in particular if A is symmetric. The uncontrolled state trajectory, also refereed to as the natural response, is then given by: where λ l are the eigenvalues of A associated to the eigenvectors v l and the coefficient vector β = (β 1 , · · · , β n ) t is given by the initial condition q 0 when solving β = P −1 q 0 , We have well-known results about controllability for discrete linear systems, they translate to CA as follows.

Proposition 2. A cellular automaton is controllable if and only if
where A describes the uncontrolled transition function and B the action of the control. If A is diagonalizable, A = P DP −1 , then the condition becomes rank[D n−1B , D n−2B , · · · ,B] = n whereB = P −1 B.
Example 4. We continue here example 3. In this case, we have that the state of our discrete linear system is given by: and q t+1 (c ii ) = F (q t (N 1 V (c ii ))) where F is deduced from (19). The matrix A corresponding to the natural response is given by: where A 1 is an m × m tri-diagonal matrix, and D 2 , D 3 are m × m diagonal matrices given by: where to simplify the notations we introduced a = α 0,0 , b = α 1,0 , c = α −1,0 , d = α 0,1 , e = α 0,−1 . A is given by a band matrix and the evolution function Evol s0 T of this uncontrolled system is determined by the eigenvalues of the matrix A. Table  1 contains the eigenvalues when assuming m = 5. It can be clearly seen by taking the limit of e and d when they tend to zero that the 25 eigenvalues collapse to the Table 1. Eigenvalues corresponding to Example 4 and m = 5.
five eigenvalues of the first column when e = d = 0. Notice that depending on the sign of the coefficients a, b, c, d, e the eigenvalues might be complex.
We define the origin of a CA has the configuration with zeroes in all grid elements, and we have the following proposition. Proof. This is a direct consequence from the fact that if the largest eigenvalue of a linear discrete control system is strictly less than 1, then the origin is asymptotically stable [3].
In particular for our example with m = 5 assume that the coefficients of the transition function are all positive, and such that a + 3bc + 3de + 6 √ bcde < 1, then the origin is asymptotically stable. In other words, the CA converges toward the equilibrium solution corresponding to a zero state value in each grid element. If a = 0.4, b = 0.2, c = 0.1, e = d = 0.1, then max i |λ i | = 0.8537 < 1 and the origin is asymptotically stable. Assume the initial configuration is such that every grid element has state 10. After 30 time steps every state for our CA is less than 0.1 and after 90 time steps it is less than 10 −6 . Choosing the parameters as a = 0.4, b = 0.2, c = 0.4, e = d = 0.2 we have that the largest eigenvalue is 1.23 and therefore strictly greater than 1. The origin is then unstable. It can be seen through simulations, indeed if we choose the same initial configuration as before (state value of 10 in each grid element) after 20 iteration all states are greater than 150, they are all increasing rapidly.
The following proposition highlights a surprising result stating that under some conditions a time invariant linear CA is controllable by controlling a single grid element only. The transition function is defined as follows:  i=1 s i (c) = 1 is always satisfied (provided it is satisfied by the initial conditions) and the state's components are positive provided that 1 ≥ α 1 +α 2 and 1 ≥ α 3 . It is a discrete model with linear and quadratic terms.
The uncontrolled fire spreading CA can be written as a discrete control system as follows: where the second terms represents the quadratic part of the transition rules, see equations (33-35). The matrix A is a block diagonal matrix: For the quadratic part, Q t = (Q t 1 , · · · , Q t 4m 2 ) is a 4m 2 × (4m 2 ) 2 with all the elements of the square matrices Q t i , i = 1, · · · , n being equal to zero except the i − th line, which equals q t and H = (H 1 , · · · , H 4m 2 ) is a (4m 2 ) 2 × 4m 2 matrix with each H i constant with mostly zeroes but some α 1 and α 2 based on the transition rules.
Remark 1. In our model, we consider N 1 V but generalizations are easy, for instance to the Moore neighborhood. Indeed, we could simply add a term in (33) (and then substract it in 35) of the form α M s 1,t (c ) with α 4 < α 2 to reflect that cells on fire that are in diagonal of the central one will have less impact for the central cell to get on fire than the adjacent ones. We also would need to adjust the condition on the coefficients to be 1 ≥ α 1 + α 2 + α 4 .
Below we provide some simulations for our CA of fire spreading. Using the red-green-blue (RGB) color model, we represent the amount of burning material with the red component, the amount of flammable unburnt material with the green component and the amount of burnt material with the blue component. Thus, the amount of inflammable material can be inferred from the luminosity of the color.
First, we compare the two different basic neighborhoods -the Von Neumann neighborhood and the Moore neighborhood. In Figure 6, we present two simulations, both of which start with a single location that is partially burning. As the two simulations progress, the different shapes that result from the two different neighborhoods become clear. It should be noted that the Moore neighborhood progresses much faster due to having more neighbors. Thus, the neighbor contribution is much greater. At timestep 50, the distance (as defined above) between the states of the two simulations is 6346.67. The timestep of the Moore neighborhood simulation that is closest to timestep 50 of the Von Neumann simulation is timestep 29 with a distance of 352.89.
Next, we compare different values of the α i 's and the effect it has on the rate of spreading as well as the width of the burning band (see Figure 7). For all the remaining fire simulations, we use the Von Neumann neighborhood. These two simulations rapidly diverge. By timestep 10 the distance between them is 147.54, by timestep 50 the distance is 4853.13 and by timestep 100 the distance is 15504.70.
For this example the control can be viewed as as an external input in the form of any treatment applied directly to burning fuel (wetting, smothering, or chemically quenching the fire) or by physically separating the burning from not burned fuel (done by fire engines, fire personnel and aircraft applying water or fire retardant directly to the burning fuel). It can also correspond to creating control lines: boundaries that contain no combustible material or by using fire retardants, fire-fighting foams, and superabsorbent polymer gels. These can clearly act on either of our three component of the state (independently or at the same time). The conditions for controllability are more complex than in the linear case, but since the matrix H is sparse the problem is easier. It will be studied in forthcoming work.
Assume we have the capability to build a barricade which for us translates into act on the fraction of the fuel cell unburned and turn it into a not flammable component. This means that we need to introduce a control that will grow s 0,t (c) by the amount s 1,t (c) and put the latter one to zero (for simplicity we do not consider where f 1 0 is given by (33-35), c ∈ N 1 V (c), c = c, and u 0,t (c) ∈ {0, 1}. If u 0,t (c) = 0 it corresponds to the natural evolution, however when u 0,t (c) is turned to 1 no more fire can take place in cell c which corresponds to the barricade being constructed.
We will consider two examples involving barricades. Computationally, a barricade is simply a region in which combustible material, burning material and burnt material are all zero. The first example, shown in Figure 8, we consider a fixed obstacle and we can see how the spreading fire has to divert around the obstacle. This increases the amount of time it takes to reach locations on the far side of the obstacles, but does not seem to have significant long-term impacts on the shape of the fire. Second, we consider the siutation where an obstacle (or sequence of obstacles) is introduced as the fire progresses. Notice that in Figure 9 a small region of fire appears to materialize on the far side of the barricade. This is due to the fact that, although it is not visible, by timestep 20, the fire has already reached beyond the barricade. After a short period of time, that small amount of fire has grown to a visible level. Creating barricades as the simulation progresses is one of the two types of control we consider.
Let us assume now that we can put the fire down in a region (which is equivalent to a subset of the cell in the grid), this means that no new fire will take place in where u 2,t (c) ∈ {0, 1} controls fire extinction. We illustrate this second form of control with a simulation starting from the same initial conditions. In Figure 10, we extinguish the fire in a small region at timestep 20. Notice that this has more significant long-term consequences on the shape of the spreading fire than obstacles. Also, notice the darker region which remains after the fire has passed. This artifact results from extinguishing an active fire, therefore reducing the total amount of material in the region.

Morphogenesis and tumor growth.
4.2.1. Cell's growth. In our model, the topographical organization of fractone expression reflects the shaping of the mass of cells. We hypothetize that the fractones guide the shape of multicellular organisms, therefore the cellular automaton must capture the interplay between the different types of cells and the fractones playing the role of chemical effector for the molecular mechanisms that control stem-cell fate. The CA in [2] was designed as two automata co-evolving in parallel and feeding each other with information to adjust the corresponding transition functions. One automaton, the diffusion automaton, was accounting for the evolution of the growth factors, and the second automaton, the mitosis automaton, was designed with a transition function depicting rules for cells division. In this paper, as a first step to develop control strategies to grow specific forms we neglect the growth factors and their diffusion process. Indeed, it was shown in [2] that the growth factor diffusion mainly controls the time scale of cell growth and that the fractones are dominant in determining shape. The main difficulty is to combine a fixed duration and a dynamic spatial location of the fractones during that duration that will bring an initial mass of cell to a desired one. The fact that the evolution has to occur over a prescribed time frame is key in morphogenesis. We introduce the state of a grid cell as a single value between 0 and 1. If 0 ≤ s(c) < 1 it is considered that no biological cell is present in the corresponding unit c of the grid, and s(c) = 1 means the existence of a biological cell. Intuitively, when s(c) ∈ (0, 1) the state expresses the probability that a biological cell will be present in the next time step. The transition rule for the uncontrolled process is for c ∈ G: where α t (c, c i ) is some constant depending on the organism under consideration. It captures the "natural" growth of a cellular mass without the action of fractones.
Here the coefficients depend on time, indeed once there is a biological cell, i.e. s t (c) = 1 the state of the cell does not evolve, unless there is an external event, which is expressed by α t (c, c i ) = 0 for all c i ∈ N (c) in that case. Note that contrary to the other example and the fire application, we have for this example a state constraint by imposing s t+1 (c) ≤ 1.
The neighborhood can be chosen in various ways depending again on the organism under study and the stage of the cellular mass under consideration (adult or embryo). In Figure 11, we illustrate a simple example of cell growth, starting from one cell and using the standard Von Neumann neighborhood. We contrast this with 100. An uncontrolled cell growth using a neighborhood consisting of the three grid units directly above the central unit and the three side-by-side units in the row three rows below the central unit. Figure 14 in which we use a non-standard neighborhood to acheive markedly different results. Although we do not illustrate this here, we also admit the possibility that the choice of neighborhood may vary with time or even be used as a control mechanism. The coefficients α are assumed such that α 2,t (c, c i ) ≥ α 1,t (c, c i ) to express that a tumor cell will divide at a more rapid rate than normal cells. The constraints make explicit that if there is a normal biological cell in grid unit c then the cell cannot become a cancer cell and can only serve to influence its neighbours to become normal cells. Also, note that N 1 and N 2 might differ since the radius of influence of an aggressive tumor might be bigger than the one for healthy cells.
The role of fractones in tumor growth is unclear at this stage, it would be however interesting to create a virtual lab to test different hypotheses with chemical composition and distribution of fractones related to the cancer cells.

5.
Conclusion. This paper is a first approach to apply discrete control theory to analyze the behavior of controlled cellular automata. The framework is designed to be flexible enough to allow for heterogeneity regarding the grid cell's states, indeed as it can be observed in the fire application for instance some areas might be set as completely inflammable while others are flammable or possibly at least partially. The heterogeneity distribution can also evolve with time, which is illustrated with the possibility to introduce obstacles throughout the evolution of the dynamical system under study. In the tumor growth application, different types of cells can be introduced to mimic organs various sensitivity to specific type of cancer.
We presented how different control strategies affect the evolution of the system, and the next step will be to use more in depth results from discrete control systems to come-up with systematic ways to design control strategies to reach a final goal. It could be stabilize the size of the rumor for instance or determine the distribution of fractones to grow a specific form during fixed time period. In addition, in forthcoming work optimal control of discrete systems will be explorer to not only reach 22 ACHILLES BEROS, MONIQUE CHYBA AND OLEKSANDR MARKOVICHENKO a prescribed goal but in the most efficient way given a cost. Based on this first approach, the authors are optimistic regarding the application of discrete control methods to this class of controlled cellular automata.