Chiral selection in wrapping, crossover, and braiding of DNA mediated by asymmetric bend-writhe elasticity

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 extra-cellular 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.

Definition 2.1. A cellular automaton is a quadruple $ CA = (S, d, N, f) $ where:

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\subseteq \mathbb 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\rightarrow \mathscr P(G) $. The most common ones are the Von Neumann neighborhood and the Moore neighborhood. Given a cell $ c\in 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\rightarrow 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.

Definition 2.2. An extended Moore neighborhood of radius $ r $, denoted $ N_M^r $, is a $ 2r+1\times 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, \cdots, k $ be the component $ p $ of the state's string. The distance between $ s_t^1 $ and $ s_t^2 $ is given by:

where $ P(c) $ is the cardinality of the set of neighbors of $ c $ over the three neighborhood $ N_M^r(c) $. It will be typically $ 35 = \sum_{r = 0}^2 (2r+1)^2 $ unless it is close to the boundary of the grid $ G $ and will therefore be smaller.

Lemma 2.4. The function defined by $ d(., .) $ on the set of configurations of a given CA is a distance.

Proof. By definition $ d $ is non-negative, and is symmetric: $ d(s_t^1, s_t^2) = d(s_t^2, s_t^1) $. It is also obvious that if $ s_t^1 = s_t^2 $ then $ d(s_t^1, s_t^2) = 0 $. Assume now that $ d(s_t^1, s_t^2) = 0 $, then $ n_{p, t, 0}^1(c) = n_{p, t, 0}^2(c) $ for all $ c $ and $ p $ which implies that $ s_t^1 = s_t^2 $. The triangle inequality $ d(s_t^1, s_t^3) \leq d(s_t^1, s_t^2)+d(s_t^2, s_t^3) $ follows directly from:

Example 1. To illustrate our metric, we use a simple example associated to a $ m\times m $ grid with single binary state $ s(c)\in \{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\in\{1, \cdots, m\} $ which means that all cells are white. It can be shown that $ d(s^0, s^2)\approx\dfrac{1}{2} $ and therefore $ d(s^0, s^1)\approx\dfrac{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)\in 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\in $ 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 $ \sim $ as follows. We say that $ s^1-s^0\sim s^2-s^0 $ if for each value associated to a cell in $ s^1-s^0 $, there exists a cell $ \tilde c\in s^2-s^0 $ with the same value and vice-versa. In other words, $ s^2-s^0 $ can be seen as a reshuffling of the cells in $ s^1-s^0 $ and their associated values.

Lemma 2.6. Let $ s^0 $, $ s^1 $ and $ s^2 $ be three configurations such that the difference $ s^0-s^1\sim s^0-s^2 $. Moreover, assume that for $ c\in s^2-s^0 $ its Moore neighborhood of maximum range considered in the calculation of the distance does not contain more than one ell with a nonzero value. Then $ d(s^0, s^1)\leq d(s^0, s^2) $. The equality holds for instance when all cell's values in $ s^0-s^1 $ (and so in $ s^0-s^2 $) are not negative (or not positive).

2.2.   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:

where $ F(s_t)(c) = f(s_t(N(c))) $.

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}\leq s_p(c)\leq c_{p, \max} $, $ c_{p, \min}\leq 0\leq c_{p, \max} $ for all $ c\in G $ and each component of the state string $ p = 1, \cdots, k $. In this case, the transition between two states is expressed as:

Example 2. The linear case corresponds to each component of the new state of cell $ c $ at time $ t+1 $ being obtained as a linear combination of the states of $ c $ and its neighbors:

where $ \alpha_{i, c'}^p:I\mapsto \mathbb 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) $.

3.   Evolution and control problem

3.1.   Evolution

We have to define the evolution (trajectory) for our system. Let $ I = \{0, 1, \cdots, T\} $ be a discrete time interval, and an initial configuration $ s_0:G\rightarrow S^G $. The evolution $ Evol_T^{s_0, u} $ for a given control $ u $ is a sequence of configuration $ s_t, t = 0\cdots, T $ given by the global dynamics:

where $ F^t $ represents $ t $ composition of the map $ F $.

Example 3. For instance, let us consider an $ m\times m $ grid and $ N(c) = N_{V}^1(c) $ to be the von Neumann neighborhood of range 1 which is composed of the central cell $ c $ and the 4 adjacent cells. If we represent a cell in $ G $ by an ordered pair of integers: $ c_{ij} $ where $ (i, j)\in \{1, \cdots, m\}\times \{1, \cdots, m\} $, then by definition we have $ N_V^1(c_{ij}) = \{c_{i+r, j+s}, |\; c_{i+r, j+s}\in G\;{\rm and}\;r, s\in \{0, -1, +1\}\} $. 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: $ \alpha_{00} = 0, \alpha_{10} = 0.24, \alpha_{-10} = 0.24, \alpha_{01} = 0.24 $ and $ \alpha_{0-1} = 0.24 $. We chose the initial configuration with ones in the diagonal and zeroes elsewhere: $ s_0(c_{ij}) = \left\{

$$\begin{array}{cc} 1\;{\rm if}\;i = j\\ 0\;{\rm if}\;i\neq j \end{array}$$ \right. $. 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)\in \{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. Figure 4 provides an example of controlled evolution with $ u(t) =

$$\begin{cases} 0 & \mbox{if $0 \leq t < 20$ , $0 \leq t < 50$ or $>t \geq 100$ }\\ 1 & \mbox{otherwise} \end{cases}$$ $ and the coefficients $ \beta_{00} = 0, \beta_{10} = -0.34, \beta_{-10} = -0.24, \beta_{01} = -1.24 $ and $ \beta_{0-1} = -0.24 $.

Section 4 provides additional examples with fire spreading and biological cell's growth.

3.2.   Controllability

Definition 3.1. A cellular automaton is said to be controllable if for any initial $ s(0) = s_0 $ and final $ s(T) = s_T $ configurations, $ T<\infty $, there exists a control $ u $ such that $ s_T = Evol_T^{s_0, u} $.

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 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\times l $ matrix and $ u_t $ is a $ l\times 1 $ vector. The matrix $ B_t $ indicates how the control is acting on the $ k $ components of the state of cell $ c $ ($ B_t\equiv 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 $ \mathbb 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\times m $, our state vector $ q $ belongs to a subset $ ( \mathbb R ^m\times \mathbb R ^m)^k $ which is determined by the state constraints. We introduce $ n = m^2k $ 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\equiv 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\times n $ matrix, $ B_t $ is a $ n\times 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 $ \phi(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 = PDP^{-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 $ \lambda_l $ are the eigenvalues of $ A $ associated to the eigenvectors $ v_l $ and the coefficient vector $ \beta = (\beta_1, \cdots, \beta_{n})^t $ is given by the initial condition $ q_0 $ when solving $ \beta = P^{-1}q_0 $, $ P = (v_1, \cdots, v_{n})^t $.

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 = PDP^{-1} $, then the condition becomes $ {\rm rank}[D^{n-1}\tilde B, D^{n-2}\tilde B, \cdots, \tilde B] = n $ where $ \tilde B = 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_V^1(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\times m $ tri-diagonal matrix, and $ D_2, D_3 $ are $ m\times m $ diagonal matrices given by:

where to simplify the notations we introduced $ a = \alpha_{0, 0}, b = \alpha_{1, 0}, c = \alpha_{-1, 0}, d = \alpha_{0, 1}, e = \alpha_{0, -1} $. $ A $ is given by a band matrix and the evolution function $ Evol_T^{s_0} $ 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 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.

Proposition 3. The origin of the CA is asymptotically stable for the uncontrolled system if and only if $ |\lambda_i|<1 $ for all $ i = 1, \cdots, n $ where $ \lambda_i $ are the eigenvalues of $ A $. It is unstable if there exists an eigenvalue $ |\lambda_i|>1 $.

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+\sqrt{3bc+3de+ 6\sqrt{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} |\lambda_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.

Proposition 4. Assume $ e = d $ and all eigenvalues of $ A $ being distinct. We also assume that $ P^{-1} $ has at least a column with non zero values. Then, the CA is controllable with a scalar control on only one of the cell's grid.

Proof. Since $ e = d $, $ A $ is diagonalizable and we assume $ \lambda_i\neq \lambda_j $ for $ i\neq j $ where $ \lambda_i $ are the eigenvalues. The rank condition $ {\rm rank}[D^{n-1}\tilde B, D^{n-2}\tilde B, \cdots, \tilde B] = n $ is equivalent to the Vandermonde matrix $ V = \left (

$$\begin{array}{ccccc} 1 &\lambda_1 & \cdots & \lambda_1^{n-1} \\ 1 &\lambda_2 & \cdots & \lambda_2^{n-1} \\ \vdots & \vdots &\vdots & \vdots \\ 1 & \lambda_n & \cdots & \lambda_n^{n-1} \end{array}$$ \right ), $ being invertible. Since det$ V = \Pi_{1\leq i<j\leq n}(\lambda_j-\lambda_i) $ it is non zero and controllability is guaranteed for any $ n\times 1 $ vector $ \tilde B $ with nonzero components. By construction $ \tilde B = P^{-1}B $ where $ P $ is the matrix of eigenvectors. If $ P^{-1} $ has a column non zero values, let us say column $ r $, then we simply set the $ r $-component of B to 1 and $ \tilde B $ satisfies the requirement.

4.   Applications

4.1.   Fire spreading

The CA model presented here is a generalization of typical cellular automaton applied to fire spreading where each cell is in only one of the following state: not flammable; the fuel cell unburned; burning cell; burned cell), and the rules are based on the neighborhood (Moore or Von Neumann) states and probabilities. Wind or external events can be incorporated in the transition rules through the use of control functions. In our model a cell can be in more than one stage with a fraction of it that can be burning and the rest of the cell already burned, moreover the external parameters can be controlled and their impact understood.

The state of a cell $ c $ is given by four values: $ s_0(c) $ is the not flammable part of the cell, $ s_1(c) $ represents the fraction of the cell that is not touched by the fire (not burning nor burned), $ s_2(c) $ is the fraction of the cell that is burning and $ s_3(c) $ the fraction that is burned. Clearly we have that $ \sum_{i = 0}^3 s_i(c) = 1 $ since it represents 100% of the cell, $ s_i(c)\geq 0 $ for all $ i $ and the not flammable part of the cell is constant throughout the evolution of the fire unless there is an external action. Therefore, we have $ S\subseteq [0, 1]^4 $.

The transition function is defined as follows:

where $ \alpha_1 s_{1, t}(c)s_{2, t}(c) $ represents the fraction of untouched cell that starts burning due to existing fire in the central cell, $ \alpha_2 \sum\limits_{c'\in N_V^1(c), c'\neq c} s_{2, t}(c')s_{1, t}(c) $ is the fraction of untouched cell that starts burning due to fire in neighbouring cells and $ \alpha_3 s_{2, t}(c) $ depicts the quantify of burning land that is now completely burned. By definition the state constraint $ \sum_{i = 1}^4 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\geq \alpha_1+\alpha_2 $ and $ 1\geq\alpha_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:

which is invertible with determinant $ (1-\alpha_3)^{m^2} $ and $ A^{-1} $ is block diagonal with $ A_1^{-1} = \frac{1}{1-\alpha_3}\left (

$$\begin{array}{cccc} 1-\alpha_3 & 0 & 0 & 0\\ 0 & 1-\alpha_3 & 0 & 0 \\ 0 & 0 & 1 & 0\\ 0 & 0 & -\alpha_3 &1-\alpha_3 \end{array}$$ \right ) $. For the quadratic part, $ Q^t = (Q_1^t, \cdots, Q_{4m^2}^t) $ is a $ 4m^2\times (4m^2)^2 $ with all the elements of the square matrices $ Q_i^t $, $ i = 1, \cdots, n $ being equal to zero except the $ i -th $ line, which equals $ q^t $ and $ H = (H_1, \cdots, H_{4m^2}) $ is a $ (4m^2)^2\times 4m^2 $ matrix with each $ H_i $ constant with mostly zeroes but some $ \alpha_1 $ and $ \alpha_2 $ based on the transition rules.

Remark 1. In our model, we consider $ N_V^1 $ 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 $ \alpha_M \sum\limits_{c''\in N_M^1(c)\backslash N_V^1(c)} s_{1, t}(c') $ with $ \alpha_4<\alpha_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\geq \alpha_1+\alpha_2+\alpha_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 $ \alpha_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 the case when only a fraction of the fuel cell unburned is turned into not flammable since it is a straightforward extension but makes the equations unnecessarily more complicated). This means that for the cells that can be affected by a barricade, the controlled transfer function is now:

where $ f_0^1 $ is given by (33-35), $ c'\in N_V^1(c), c'\neq c $, and $ u_{0, t}(c)\in \{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 such cell at that moment and the fire will go extinct. However, it does not convert the "fuel" fraction of the cells into not flammable. The control would then be introduced as follow:

where $ u_{2, t}(c)\in \{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.

4.2.   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 \leq 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)\in (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\in G $:

where $ \alpha_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 $ \alpha_t(c, c_i) = 0 $ for all $ c_i\in 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)\leq 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 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 natural growth without the introduction of fractones will provide a growth directed by the constants $ \alpha_t(c, c_i) $, in particular if we start with one cell and we assume the coefficients $ \alpha_t(c, c_i) $ do not depend on $ c_i $ but are all the same for every cell in the neighborhood we will grow a isotropic mass (see Figure 11).

To break the symmetry and create more sophisticated masses of cells we must activate the controls (playing solely with the notion of neighborhood is not fine enough and does not reflect our biological hypothesis). The controls are: the fractones that can accelerate or decelerate growth; and spontaneous death for the cells.

Definition 4.1. If a fractone is associated to a cell, then its corresponding coefficient in the transition function is altered either by an increasing factor if the fractone is inducing accelerated growth, in which case we call it a positive fractone, or by a decreasing factor if the fractone is decelerating growth which is then called a negative fractone.

More precisely, we have the following. Assume $ c_j $ such that $ s_t(c_j) = 1 $, i.e. there is a biological cell in that grid unit, and that we have a positive fractone associated to this biological cell. We introduce $ u^+_t(c_j) $ a positive real number and the coefficients $ \alpha_t(c, c_j) $ are replaced by $ \alpha_t(c, c_j)+u^+_t(c_j) $ in all transition rules that involve $ c_j $. In particular, if $ c_j $ is in the neighborhood of $ c $, we have:

The value $ u^+_t(c_i) $ depends on the chemical composition of that specific fractone. The transition rule is modified the same way for negative fractones with the constraint that $ \alpha_t(c, c_j)-u^-_t(c_j)\geq 0 $, $ u^-_t(c_j) $ is a positive real number.

The present model can be described in terms of a discrete control system as follows. Assume the possible chemical composition of positive fractones is captured by $ l $ positive real numbers (i.e. there is a finite number of speed at which fractones can accelerate growth), and similarly for negative fractones we assume their possible chemical composition is given by a finite vector of length $ s $: $ \beta^+\in \mathbb R ^l, \beta^-\in \mathbb R ^s $. Assume the dimension of the state vector $ q $ to be $ n $, we then have:

where $ A_t $ is an $ n\times n $ matrix that represents the natural growth, and $ U^+ $ (resp. $ U^- $) is a $ n\times l $ (resp. $ n\times s $) control matrix whose entries are binary, 0 or 1, with 1 representing the existence of a positive (resp. negative) fractone. Notice that the control part does not include the state since fractones are associated to existing cells and therefore only to $ c_j $ such that $ s_t(c_j) = 1 $. The value of $ u^+_t(c_j)) $ ($ u^-_t(c_j) $) is deduced from $ U^+_t\beta^+ $ ($ U^-_t \beta^- $).

In Figure 13, we illustrate the affect that the fractones can have on the shape of the results. One line of fractones increases the rate of growth and results in the oblong shape of the cell mass. Below that is a line of fractones that stop cell growth. Notice that these hardly have any impact on the shape, except that they stop cell growth at their locations. Above the first line is a row of blocks of fractones that stop cell growth. These have some impact on the shape of the cell mass simply because they act as an barricade, much like the barricades in the fire simulations.

Definition 4.2. Death of cell can occur in a spontaneous way. If a biological cell exists in grid unit $ c $, i.e. $ s_t(c) = 1 $, the death of this specific cell means that its state reset to the value zero: $ s_{t+1}(c) = 0 $.

By definition the control corresponding to cell's death is a vector $ u^d_t\in \mathbb R ^n $ whose entries are zero for state components $ q_t^i $ that are not equal to 1 or correspond to an existing biological cell not dying, and are 1 for grid units with a biological cell undergoing spontaneous death. The system becomes:

4.2.2.   Tumor growth

A single grid cell of our two dimensional lattices will be occupied by either an empty space, a cancer cell or a normal cell (not cancerous). The state for each grid cell is characterized by an ordered pair $ (s_1, s_2)\in [0, 1]^2 $ to mimic our model on cell's growth 4.2.1 where $ s_1 $ accounts for the probability to get a normal cell and $ s_2 $ for the probability to get a cancer cell. If $ s_{1, t} = 1 $ there is a biological cell in that grid unit and if $ s_{2, t} = 1 $ there is a tumor cell.

The model is similar to that of cell growth, but a state is now a pair and the transition rules for the natural growth are as follows:

with the additional state constraints:

The coefficients $ \alpha $ are assumed such that $ \alpha_{2, t}(c, c_i)\geq \alpha_{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 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.