Co-evolving cellular automata for morphogenesis

Morphogenesis, the shaping of an organism, is a complex biological process accomplished through an well organized interplay between growth, differentiation and cell movement.It is still today one of the major outstanding problems in the biological sciences. Pattern formation has been well-addressed in the literature with the development of many mathematical models including the famous reaction-diffusion ones. We here take a different approach, introducing a controlled cellular automaton in order to model the signal molecules, known as growth factors, that convey information from one cell to another during an organism's development and help maintain the viability of the adult. This control represents extracellular structures that have been associated with the regulation of stem cell proliferation and are called fractones. In this paper we introduce two co-evolving automata, one describing the perturbed diffusion of growth factors and one accounting for the rules of basic cellular functions (proliferation, differentiation, migration and apoptosis). Fractones are introduced as an external input to control the shaping of multi-cellular organisms; we analyze their influence on the emerging shape. We illustrate our theory with 2 and 3 dimensional simulations. This work presents the foundation upon which to develop cellular automata as a tool to simulate the morphodynamics in embryonic development.


1.
Introduction. An important question in biology is how the developing embryo is organized on a cellular level to grow into a specific morphology (shape). Growth occurs through mitosis of individual cells, the splitting of one mother cell into two new daughter cells, and it is only through repeated acts of mitosis over time that a cell mass can form. While some influencing factors have been identified (forces between cells [5,19] and growth factors detected by the cells through receptors on the cell surface [24]), it is still uncertain how the embryo organizes and directs cellular growth. In his pioneering work on mathematical modeling of morphogenesis [22], Alan Turing introduced his now well-known intercellular reaction-diffusion process. His model is based on a system of chemicals diffusing in the extra-cellular space and reacting with each other. In the model presented here we consider specialized extra-cellular matrix structures, which have been named fractones [15], as a refined control mechanism orchestrating the coordinated diffusion of the growth factors through the extra-cellular space.
We focus on four principal cell functions which influence the development and viability of cellular masses: Proliferation; Differentiation; Migration and Apoptosis. We base our model on the fact that these four basic cell functions are regulated by growth factors which translate signals into cell functions. Some stimulate cellular division (these can be specific to a particular cell type or not) and some stimulate cell and tissue function through influencing cell differentiation. For instance, it has been shown that the fibroblast growth factor receptor 2-IIIb and the vascular endothelial growth factor receptor 2 can promote cell migration and proliferation [13]. Peptide growth factors (and Bcl-2 gene), first characterized as stimulators of cell proliferation, are now also known to inhibit death [6,7].
Existence of a governing mechanism is key to providing the specificity and precision that the diffusion of chemical signals needs for proper development. Recently discovered in the neurogenic zone of the adult mammalian brain [15,17], fractones are microscopic extracellular structures that have been associated with the regulation of stem cell proliferation [9,12,16]. See Fig. 1 for an illustration of fractones around the aqueduct obtained using immunofluorescence labeling of fractones and imaging with confocal laser scanning microscopy. In particular, it has been shown that fractones can promote cell proliferation [9] and can also inhibit cell proliferation [8] in the adult brain. More precisely, since the basal laminae is thought to contribute to morphogenesis, it was hypothesized that fractones were directly influencing the cells with which they formed connections. This hypothesis was supported by experiments detecting the existence of heparan sulfate proteoglycans (HSPG) in fractones which capture growth factors (specifically FGF-2) for presentation to cells. The fractones are thought to capture and store these growth factors until this concentration reaches a prescribed threshold to initiate binding with cell-surface receptors. Importantly, not all fractones were found to have HSPG, suggesting that there may be specialization among fractones as to the types of growth factors they will bind [12]. It is thought 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 goal of this paper is to use cellular automata to model the hypotheses described above. Cellular automata (CA) are well-suited to mathematically capture the essence and complexity of biological systems [2,10]. They provide an efficient method for producing sophisticated self-organized structures representing systems of simple agents which, through their interplay, exhibit complex intelligent behavior [11,23,25]. Our paper is a continuation of an initial approach presented in [4] where co-evolving automata modelling diffusion and mitosis were introduced. In this work, we expand the model to incorporate how multiple growth factors, possibly triggering antagonistic actions, coordinate to determine stem cell fate by targeting selected stem cell populations [15]. Moreover, the CA has been expanded to three dimensions and we have developed in parallel a 3D-vizualisation tool used for the illustrations in this paper. The simulations presented in this paper are done on a simplified version of the introduced model, indeed the goal was to illustrate some basic properties of the evolution of the automaton to single-out the role of the controls. Forthcoming work will address numerical simulations that will take into account multiple cell types and fractones with different chemical feature to incorporate differentiation, migration and apoptosis. The outline is as follows. In Section Figure 1. Fractones are visible as the small puncta arranged around the aqueduct. The larger elongated shapes surrounding aqueduct are blood vessels. Image of the cerebral aqueduct in an adult mouse taken by 20x PlanApo objective lens. The aqueduct connects the third and fourth ventricles together. The image is from [1], see [18] for information about the experimental technique. It was obtained in Dr. Mercier's laboratory at the University of Hawaii.
2 we introduce the two co-evolving automaton with their rules and describe their interaction. We start with a two-dimensional version for simplicity and then expand to a three-dimensional one. Section 3 introduces two metrics as well as the control mechanisms and problem associated with our mathematical and biological process. Section 4 illustrates the theory and the evolution of the co-evolving automata and Section 5 discusses briefly future work as well as possible validation for our model.

2.
Cell-fractone automaton. It is well known that Cellular Automata (CA) are well-suited to produce sophisticated self-organized structures using simple agents which, through their interplay, exhibit complex intelligent behavior [11,23,25]. A CA is described by a finite collection of states for each unit in the grid on which it is defined (we restrict the word "cell" for biological cells in this article) and a function (or transition rules) that describes how the state of a unit should change at the next time step based on its current state and that of some neighboring states. Transition happens at each unit simultaneously in discrete time steps. Transition for a cell is here determined by the biological rule that when a fractone has accumulated enough growth factor corresponding to an action then the corresponding action takes place, it is therefore not random process as it can bee seen in many application (see [3] for instance). This is justified by the hypothesis that morphogenesis is not random, it is a well controlled mechanism with actions taking place when specific conditions are met and not under probabilistic rules.
Our primary goal is to determine how, in our model, the topographical organization of fractone expression reflects the shaping of the developing brain. More precisely, 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 phenomenon is one of a perturbed diffusion for the growth factors due to the presence of fractones. These act like sinks and depending on their chemical composition capture specific growth factors to either promote or inhibit cell functions. The cell functions that we consider in our model are: Proliferation: Cell proliferation is the process by which cells reproduce themselves by growing and then dividing into two identical copies. An important component of the cell cycle in cell proliferation is mitosis, when replicated chromosomes are separated into two new nuclei followed by the cell's division. Differentiation: Differentiation is the process by which a generic embryonic cell becomes a more specialized cell type. During the development of a multicellular organism a cell might go through differentiation multiple times. Migration: Cell motility is a fundamental cellular function and is key in physiological processes including development and tissue repair of multicellular organisms. Migration also plays an important role in some pathological processes like cancer metastasis. Apoptosis: A multicellular organism is a highly organized community and its number of cells is tightly regulated not only through proliferation but also by controlling the rate of cell death. Apoptosis is a form of suicide that takes place by activating an intracellular death program.
Our approach is to define two co-evolving automata: the cell-functions automaton CA c defined by the rules according to the four basic functions introduced above; and the diffusion automaton CA d that captures the perturbed diffusion of the growth factors. Below we describe first the co-evolving automata in two dimensions and then expand to three dimensions.
2.1. Two-dimensional automaton. The automata are defined on superposed grids with the diffusion one being more refined. Indeed, for the diffusion automaton a unit corresponds to a fractone. In order to respect the biological scale, this implies that a cell is a 9 × 9 region in that discretization. This is based on measurements of developing rat embryos which show the neuron soma size to be approximately 72 − 81 µm 2 in size [14] and the averaged 9 µm cell diameters of yeast cells [26]. Fractones are therefore assumed to have a diameter of 1 µm, which is on the smaller end of the 1-4 µm diameter found experimentally. Assuming a channel of 4 units in between the cells, we have that the discretization of the mitosis automaton is based on units that cover 13 × 13 units in the diffusion discretization. See Figure 2 for an illustration of the grids and a sample cell-fractone configuration.
For our model, we consider mitosis to be an instantaneous action in the sense that a cell immediately divides into two new cells. We account for the time mitosis takes by enforcing a delay before the new cells can undergo mitosis. We also neglect the forces in play during migration (this process occurs with a low Reynolds number) and assume instantaneous migration with the initial position becoming empty and the final position having a new cell. Taking cell deformation and forces into account would require an automaton with a non uniform grid and it would add too much complexity at first. Differentiation and apoptosis are also considered instantaneous events. Next we define the rules of the co-evolving automata.
2.1.1. Diffusion automaton. We introduce G d the diffusion grid, and a function N d which, for each u ∈ G d , returns the neighborhood of u. We assume here that the neighborhood consists of the eight units surrounding a central unit, i.e. for u ∈ G d : where N, E, S, W represent the four cardinal directions.
Definition 2.1. The state of a unit in the diffusion grid takes its value in {0, 1, 2} where 0 corresponds to an empty unit, 1 to a fractone and 2 to a unit in a biological cell. We will denote by s d : G d → {0, 1, 2} the diffusion state, and we introduce Note that for the diffusion automaton the type of cell is not important, we can treat them all the same since the role is solely to be an obstacle for the diffusion.
As explained in the introduction, growth factors are key elements to the process and the diffusion automaton will capture their perturbed diffusion. Each growth factor has the capability to bind with receptors associated with all four basic functions, however, as mentioned in Section 1, a specific growth factor is specialized to promote/inhibit some functions in priority. For instance FGF family members bind heparin and possess broad mitogenic and angiogenic activities. To capture this feature in our model we associate to each growth factor a 4-tuple whose components are the required concentration thresholds for binding with the receptors for each cell function. Using this system, we would describe a highly effective signaling molecule by having a low threshold, implying that very little of the growth factor is needed to trigger a cell action.
1. Each unit in the grid is assigned a list of non negative real numbers representing the concentration of the various growth factors in that specific unit. It is represented by a function g : G d → R m ≥0 , with g(u) = (g 1 (u), · · · , g m (u)) t where g i (u) is the concentration of the i th growth factor in unit u and m the number of growth factors under considerations. We refer to g as the growth factor concentration. We impose in the sequel that g i (u) = ∞ if s d (u) = 2 to reflect that biological cells are regarded as obstacles for the diffusion (by construction each obstacle is a 9 × 9 square region). 2. To each function g i , we associate a vector is the concentration threshold corresponding to proliferation (resp. differentiation, migration and apoptosis). It is assumed to be constant throughout the evolution of the automaton. Definition 2.3. We define the chemical composition of a fractone as its sensitivity to a prescribed growth factor. Let us consider a grid element u ∈ G d corresponding to a fractone (i.e., s d (u) = 1). We associate to u a list of m real positive numbers between 0 and 1, and define the fractone sensitivity f : represents a percentage describing the concentration of growth factor i captured by fractone u.
Unlike the growth factors, in our model the chemical composition of fractones (i.e. its fractone sensitivity) can possibly vary throughout the life of a fractone (both the existence and sensitivity of a fractone are determined as external inputs and not in the rules of the automaton, see Section 3).
The rules for CA d can now be expressed as follows.

Rules.
1. For a growth factor i, consider u ∈ G d . Case s d (u) = 0: If the neighborhood contains one or more fractone neighbors, i.e. ∃e ∈ N d (u) such that s d (e) = 1 (possibly more than one), then the transfer will always be to a fractone rather than any other location.
In this case, a fractone e ∈ F d in N d (u) is selected randomly according to a uniform discrete distribution and f i (e) units of growth factor are transferred from u to the fractone. Otherwise, we select the grid elements in N d (u) such that their growth factor concentration is less than in unit u, i.e. we construct the set: If L i (u) = ∅, then no transfer of growth factor i occurs. If not, one unit of growth factor i is transferred to an element of L i (u) randomly selected according to a uniform discrete distribution. Note that if s d (e) = 2, we have by definition that g i (e) = ∞ which prevent growth factors from ever being moved into the obstacle locations. Case s d (u) = 1: The grid element u is a fractone and no transfer occurs. Case s d (u) = 2: The grid element u belongs to a biological cell and no transfer occurs. 2. When a unit in the cell function automaton undergoes mitosis, a new cell is created causing a new obstacle to be created for the diffusion automaton (see Section 2.1.2 for more details). If a cell migrates, its former position becomes empty and is available again for diffusion of the growth factors while the new position is now a new obstacle as for mitosis. In case of cell apoptosis, the previously occupied space is now empty and diffusion can occur. In the case of an empty space being filled with a cell, the state of the corresponding units changes from 0 to 2 and the state of units linked to apoptosis or disappearance of cell due to migration changes from 2 to 0. Finally, the growth factors in the 9 × 9 region newly occupied by a cell (either due to proliferation or migration) are displaced to the immediate boundary of the region following a straight line along one of the cardinal or intercardinal directions (again using a uniform discrete distribution among the eight possible directions).

2.1.2.
Cell-function automaton. The cell-function automaton evolves based on information fed from the diffusion automaton. As explained earlier, the discretization of the mitosis automaton is based on units that cover 13 × 13 units in the diffusion discretization (see Figure 2), and we denote by G c the cell-function grid. In the sequel, we say that a fractone is associated with a biological cell if in the diffusion grid G d it is located in a unit that belongs to the border of the 9 × 9 cell's region. When simulating mitosis, an important question is where daughter cells will be located. A natural approach would be to consider empty spaces for the mass of cells to grow within the Von Newman neighborhood in the G c grid which would be formed of the four adjacent units around a central unit. However, this eventually limits mitosis to cells located on the border of the mass of cells under consideration only which does not reflect the biological reality. An alternative would to look at empty spaces without any limitation of distance around a central unit and permit growth in these directions using a pushing algorithm. We hypothetize that this again will eventually not reflect the biological reality. Therefore, in preparation of future work to account for different types of cells and communication channels between them, we introduce an extension of the simple von Neumann neighborhood by considering units at a Manhattan distance (see Figure 3), less than or equal to r which is a fixed number to be determined: Definition 2.4. The state of a unit in the cell-function grid takes its value in {0, 1, · · · , k} where 0 corresponds to an empty unit, 1 to a unit that contains a stem cell and numbers between 1 and k refer to different types of cells (muscle cell, blood cell, nerve cell, cardiac cell, liver cell, etc.). We denote by s c : G c → {0, 1, · · · , k} the cell-type state. We introduce M c = {u ∈ G c : s c (u) = 0} ⊂ G c the mass of cells, and M i c = {u ∈ G c : s c (u) = i} ⊂ G c , i = 0 the i-type mass of cells. To induce one of the four cell functions under consideration in this paper to occur, the hypothesis is that a fractone stores a quantity of a growth factor until it reaches the threshold concentration corresponding to this specific cell function. The likelihood for a fractone to trigger a specific event (proliferation, differentiation, migration or apoptosis) is therefore an interplay between the fractone sensitivity and the concentration thresholds ω i . We now define the notion of activity for fractones.
Definition 2.5. With respect to growth factor i and a specific function v (v ∈ {p, d, m, a}), a fractone becomes active once the threshold concentration ω v i is captured by the fractone. It is then capable of acting on the fate of the cells. If a fractone has not reached a threshold concentration for any growth factor it is said to be passive.
Note that while a fractone might have reached the required threshold concentration, it is hypothesized that it might choose to defer the time at which it will induce action. As a first approach and to simplify the complexity of the algorithm, we will assume in this paper that action follows immediately once the threshold is attained.
Definition 2.6. To each biological cell we associate its age. This is done by tracking the time steps during the evolution of the mitosis automaton. More precisely, we define the age function as: When mitosis occurs, a mother cell gives birth to two daughter cells and the counting of time steps is reset to 1. However, when differentiation takes place we do not reset the age of the cell.
Rules. At every iteration step, CA c queries CA f for the existence of an active fractone. If such fractone is identified the corresponding unit is analyzed using the following rules.
Proliferation: If u, a unit in the cell function grid, is associated with an active fractone for proliferation there are two possible outcomes. Either the growth factor promotes division or inhibits it. 1. If the growth factor promotes proliferation and at least one neighboring unit is empty, then mitosis will occur. The algorithm decides which neighboring empty unit will change state (i.e. will have a new cell) randomly using a uniform discrete distribution among the ones at the closest manhattan distance from u. If a central unit is associated with more than one positive active fractone, one is selected randomly. The selected fractone's growth factor concentration is then reset to a fraction, σ, of its previous concentration. Because multiple cells might have a shared empty neighbor, a mitosis order needs to be designed to avoid conflicts. In this paper, we use the lexicographic order within our grid. As a consequence, a unit ready for mitosis at the beginning of a timestep may not undergo mitosis if all empty neighbors are taken by other units. 2. If the growth factor is an inhibitor for proliferation it slows down the next mitosis event. More precisely, when a fractone is active for inhibition of proliferation it awaits the event to become active as well for promoting proliferation (a given fractone captures multiple growth factors) and defers the scheduled mitosis for a pre-defined number of time steps. At that stage, the inhibitor growth factor concentration is then reset to a fraction of its previous concentration. Differentiation: When a cell changes type, the cell-type state function is affected and the threshold concentration of the corresponding active fractone is reset to a fraction of its previous value for the growth factor that was used to trigger the chemical reaction. Migration: This central process requires multiple active fractones as movement of cells has to be orchestrated by more than one factor. We model that using fractones in the following way. The migrating cell will have an active migration fractone signaling a request for motility, and a cell at the desired location of migration will also have an active migration fractone acting as a "magnet" to attract the other cell (we assume that different growth factors can generate different migration actions, either motility or attraction). If more than one fractone is activated as a magnet for the migrating cell, the former will choose randomly the migration location among the possible ones. Apoptosis: If a cell is no longer needed, suicide is committed by activating an intracellular death program. If a fractone is active for apoptosis, then the cell simply disappears as well as any fractones attached to it.
If there are active fractones associated with a cell that trigger different effects, then the order in which the fractones act is randomly chosen. The order may prevent some of the fractones from acting -for example, if apoptosis occurs before differentiation.
Fractones have not been observed migrating, therefore we make the assumption that fractones stay in the same place after any cell function takes place and disappear in the case of apoptosis or migration of the associated cells. The rules for the fractones are set by the external user, see Section 3.

2.2.
Three-dimensional automaton. The 3-dimensional cell-fractone automaton is a straightforward expansion of the 2-dimensional one. The rules are almost identical modulo some adjustments listed below.
Rules. We here describe the needed modifications to go from two to three dimensions.
Grid: A unit is now a cube in the 3D-grid and a biological cell is a 9 × 9 × 9 region in the discretization of the diffusion grid. Diffusion Automaton: 1. For a grid element u, we now need to consider a 26-cell cubic neighborhood N d (u). 2. The displacement of growth factors due to the creation of a new cell after mitosis is similar with the new concentration distributed among the one unit channel around the boundary of the newly covered region by the cell.
There are now 26 possible directions versus 8 in two dimensions.

Cell-function Automaton:
The 3D neighborhood N c (u) is the straightforward extension from the two dimensional neighborhood which uses the 3D Manhattan distance instead of the 2D Manhattan distance (see Figure 4). The rest of the rules are identical. 3. Cell-growth control. In our model, the rules of the co-evolving automata do not capture the totality of the biological phenomena associated with stem cell fate. Indeed, neither the creation/disappearance of fractones nor the production of growth factor are accounted for by the rules of our model. The reason is our lack of understanding of the biological processes as well as the desire to study how morphogenesis is affected by fractone population, chemical sensitivity and localization. Introducing the fractones as an external input rather than within the rules of the automaton allows us to control the development of multi-cellular organisms as well as introduce and study possible anomalies (numbers significantly larger than what can be observed on biological maps, or increased sensitivity with respect to a particular cell-function).
In other words, we have a control system with two controls to guide the shaping of an organism. The first control is the production, in terms of source and intensity, of the growth factors. Note that these chemical substances have long been associated with morphogenesis, however, as it was mentioned in the introduction a simple diffusion cannot create the complexity of shapes observed in living organisms. The second control in our process, the fractones, precisely perturb the diffusion by acting as captors and provide a more sophisticated and refined guiding mechanisms for the morphogenetic events. To introduce a framework for our control system we need first to introduce a metric.
3.1. Morphogenetic metrics. We introduce two notions of distance between configurations of our cell-fractone automaton. A configuration C is formally defined on G c by its cell-type state s c and its age function: C = (s c , age). 1. The geometric distance between two configurations C 1 , C 2 is defined by the Hausdorff distance using the Manhattan metric d 1 on the cell-function grid between the corresponding two masses of cells M 1,c and M 2,c : Note that the distance d g solely measures how the shapes of two masses of cell differ without consideration for the type of cells, forthcoming work will consider variations of this distance that account for the type of cells present in the mass as well. To measure the biological proximity of two masses of cell we must also take into account the age of the cells.
Definition 3.2. The biological distance between two configurations is defined as: The motivation for this choice of distance is to compare the age of cells that are close physically, if two cells from the two different masses are distant on the grid their age difference will most likely not be taken into account. Remark 1. The placement of the fractones plays a critical role since it is relevant for the fate (growth or stability) of the biological organism. For instance two masses of cell with geometric distance zero will evolve completely differently if one of them is void of fractones and the other one is at the opposite associated with a large number of fractones. However, since we consider them as the control in our system they can be instantaneously altered to match a prescribed distribution which is why they do not appear in our metrics.

Control mechanisms.
Definition 3.3. An evolution interval, I = {1, · · · , T }, is an interval for some length T of timesteps under consideration for the evolution of an automaton.
3.2.1. Growth factors as controls. The growth factors concentration in units of the diffusion grid are represented by the growth factor function g. There are two possible ways this function can get altered as time evolves. First through diffusion, this is described by the rules of the diffusion automaton. In addition, we can add sources of production of growth factors. These are defined on G d as well and can be localized to accelerate actions in some regions for a given mass of cells. Experimental maps do not provide information about the concentration and production of growth factors, therefore we will consider the sources of production of growth factors as a control to study how they affect the evolution of the mass of cell.
Definition 3.4. We define the growth factor control as a function w growth defined on I such that w growth (k) : G d → R m ≥0 . During an evolution interval of the cell-fractone automaton, if there exists k ∈ I such that w growth (k) is not the zero function, then the growth factor concentration is instantaneously affected at time step k and we have g(k) → g(k) + w growth (k) before the rules of the automaton are applied.

Fractones as controls.
At every time step, the position of the fractones needs to be provided as an external input. The biological mechanism creating and destroying the fractones is not well-know and therefore cannot be included as a rule. Moreover we would like to analyze the dependance of morphogenesis with respect to the spatial distribution of fractones throughout the evolution.
Definition 3.5. We define the fractone control as a function w f ractone defined on I such that w f ractone (k) is a set of pairs. The first component in each pair is taken from a subset of the units in G d such that their state is 1 and the second component is a fractone sensitivity. This determines the set of fractones F d (k) at each iteration k during the evolution.

Trajectory.
We can now define a trajectory for the cell fractone automaton. Definition 3.6. A cell-fractone trajectory defined on an evolution interval I = [0, T ] is a sequence τ = (τ (k)) T k=1 such that: where s k c , age k , g k are respectively the cell-function state, the cells' age and the growth factor concentration at time step k obtained from applying the rules of the cell-fractone automaton with corresponding controls w growth and w f ractone and issued from an initial value τ (1).
Note that the diffusion state follows from the cell-function state and the control function placing the fractones.
Remark 2. The controls cannot be pre-determined to produce a specific trajectory (i.e., cannot be an open-loop controller) since the controls depend on the evolution of the diffusion space and the evolution of the cell mass, both of which contain some degree of randomness. They are also not closed-loop controllers since they are not completely determined by state of the automaton as they include an external controller. The more inputs we obtain from the experimental biological maps the more closed-loop the controllers will become.

Control problem.
A natural control problem to consider is the following one. Given an initial (respectively final) growth factor distribution and mass of cells (including their type and age) τ initial (respectively τ f inal ), where T is fixed, can we find controls w growth and w f ractone such that its corresponding cell-fractone trajectory (τ (k)) T k=1 satisfies τ (1) = τ initial and τ (T ) = τ f inal . This is restrictive since we are asking to match not only the cell-function states but the growth factors concentrations as well and the latest cannot be obtained by experimental maps. Moreover the concentration of growth factors can be easily and quickly altered through the growth factor control by introducing new sources. For this reason, we relax controllability on the growth factors concentrations and require that we match the cell-function state and age of the cells only. We remind the reader that the localization and sensitivity of the fractones belong to the control and can be modified at every time step.
1. The cell-fractone automaton is said to be controllable in time T < ∞ from C 1 to C 2 if there exist controls w growth , w f ractone defined on I T , as well as an initial growth factor concentration g initial such that the corresponding cell-fractone trajectory τ satisfies where M i 1,c (resp. M i 2,c ) is the mass of cells of type i in C 1 (resp. C 2 ), age 1 (resp. age 2 ) is the age function associated with configuration C 1 (resp. C 2 ).
It is said to be weakly controllable in time in other words, we remove the age requirement. 2. The cell-fractone is biologically controllable in time It is said to be weakly biologically controllable in time T < ∞ if given ε > 0 Remark 3. The introduction of the notion of biologically controllable comes from the fact that there is randomness associated with our rules which makes controllability virtually impossible. This is aligned with the fact that biological shapes with a large number of cells evolved from a single initial cells and under the same rules might present slight variations.
There is a lot of flexibility with the fractone automaton, it is actually quite straightforward to prove weak controllability using fractones and production of growth factors in a suitable way. Indeed, starting from an initial mass of cells and using a single fractone we grow the shape by one cell at a time staying as close as possible in the boundary of to the final mass of cells and not exceeding the number of cells in the final mass. Then we modify fractones and sources of growth factors to induce migration of cells (possibly apoptosis in case we did grow too many cells) and differentiation. We can that way create exactly any final biological shape. Controllability is not as clear since the age of the cells need to match which imposes a much more structured way to shape the growing organism. This said, the main difficulty in the controllability problem for biological systems (whether it is any definition introduced above) is the constraint imposed on the duration of the process. Morphogenesis is a very regulated process that takes place over well defined time constraints, and this part is extremely challenging to address. Control theory on cellular automata is under-developed and it is done mostly in an exploratory way based on the specific automaton under study.

4.
Simulations. In this section, we provide 2D and 3D simulations for our fractone automaton to demonstrate the evolution of the shaping organism under our rules. The simulations and data provided in this paper are limited to a unique growth factor and one cell function, namely proliferation, to explicit the direct correlation between the two controls and the evolution of the organism. The growth factor is assumed to promote division and each fractone is set to have sensitivity 1 to this growth factor. 4.1. The role of w growth (k). The following data illustrate the impact of the production of growth factor on the time evolution of cellular activities during morphogenesis. We produced three simulations all of which began with a single cell with one fractone attached to the cell. The fractone position did not change throughout the simulation. The difference between the three simulations was the intial distribution of the growth factor. The initial distributions were as follows.
1. A uniform distribution of growth factor between 1 and 10 at each unit in the diffusion grid. 2. A plane of very high concentration near the cell on the opposite side of the cell from the fractone (the cell is, therefore, acting as a barrier to growth factor reaching the fractone). 4.2. 2D simulations of starfish-like shape. In Figure 5 the co-evolving automata produce a shape with a center and legs starting from four cells positioned into a square. This simulation clearly illustrates the impact of w f ractone on the evolution of the shape. The initial configuration is four adjacent cells, each having a single attached fractone. We control the shape by associating a direction to each fractone, meaning that the fractone promotes growth in that direction; at each stage that a new cell was created that was further in the associated direction, the original fractone was destroyed and a new fractone was created attached to the new cell. At the start of the simulation, each unit in the diffusion grid was assigned a concentration of growth factor between 0 and 20 according to a uniform discrete distribution. The six snapshots shown are at timesteps 1, 35, 100, 130, 160 and 198. In the final timestep, the oldest cell is 195 timesteps old and the average cell age is 96 timesteps.
4.3. 3D uniform growth. Figure 6 shows snapshots of a cell mass that developed from a unique cell and using a single fractone throughout the simulation. The fractone is associated with the original cell and stays at the same position there the entire time. To ensure cell mitosis within a reasonable time, the simulation is initialized with a high concentration of growth factor at every point of a horizontal plane just below the starting cell (in red in the figure).
4.4. 3D elongated growth. With Figure 7 we provide a 3D illustration of how the control w f ractone can be used to direct the growth in a particular direction, thereby creating an oblong cell mass. For this simulation, we initialize the growth factor with a uniform discrete distribution at every unit of the diffusion grid. During the whole simulation there is at most one fractone, and that fractone dies and a new one is created each time there is a new cell further in a fixed direction. During the course of the simulation, the existing fractone dies and a new one is created 12 times. At the final stage, there are 41 cells and their average age is 115 timesteps.

5.
Conclusion. This paper presents a first approach using a cellular automaton based on the existence of fractones to model morphogenesis. It was done by introducing two co-evolving automata feeding information to each other. Simulations provided an illustration of the perturbed diffusion process and some of the basic rules for cellular functions. It is important for the reader to realize that while the model possess a lot of flexibility in terms of controllability, nature imposes the time frame to grow a specific shape. This prevents for instance to grow a shape one cell at a time. The control problem is therefore quite complex. While the current simulations only include mitosis, forthcoming work will include differentiation to account for multi-cellular organisms as well as the essential functions that are migration Figure 6. 3D simulation of a uniform mass starting from a unique cell and using a unique fractone throughout the simulation. The cell color indicates age -a red cell is the product of a recent mitosis, a yellow cell is old enough to undergo mitosis. The color of the diffusion units indicate concentration with blue representing low concentration and red high concentration. Note that for 3D simulations, we have hidden the diffusion concentrations except close to each of the cells in order to make the shape of the cell mass visible. and apoptosis. A major focus will be to understand how to optimize the cellular automaton to create organisms with a much larger number of cells. Validation of the model is vital; for this we plan on running comparisons between our simulations and biological maps obtained using a confocal laser scanning microscope. The end of the paper introduced notions of distance and control problems associated with morphogensis. This is clearly highly nontrivial and the literature on systematic ways to develop control strategies on cellular automata does not exist to our knowledge (see [20,21] for some attempts) and needs to be developed. This will be the main focus of our forthcoming work, especially in the context of optimality or viability.