MORPHOGENESIS MODELIZATION OF A FRACTONE-BASED MODEL

. It has been hypothesized that the generation of new neural cells (neurogenesis) in the developing and adult brain is guided by the extracellular matrix. The extracellular matrix of the neurogenic niches features specialized structures termed fractones, which are scattered in between stem/progenitor cells and bind and activate growth factors at the surface of stem/progenitor cells to inﬂuence their proliferation. We present a mathematical control model that considers the role of fractones as captors and activators of growth factors, controlling the rate of proliferation and directing the location of the newly generated neuroepithelial cells in the forming brain. The model is a hybrid control system that incorporates both continuous and discrete dynamics. The continuous dynamics of the model features the diﬀusion of multiple growth factor concentrations through the mass of cells, with fractones acting as sinks that absorb and hold growth factor. When a suﬃcient amount has been captured, growth is assumed to occur instantaneously in the discrete dynamics of the model, causing an immediate rearrangement of cells, and potentially alter-ing the dynamics of the diﬀusion. The fractones in the model are represented by controls that allow for their dynamic placement in and removal from the evolving cell mass.


1.
Introduction. The focus of this research is to construct a mathematical model of morphogenesis, the biological process by which we take on our particular shape. Our model seeks to study how the developing embryo is organized on a cellular level in order to grow into a specific morphology. Our growth occurs through mitosis of individual cells; it is only through repeated acts of mitosis over time that the cell mass can form. This growth, however, is not random, and so must be controlled by a governing mechanism.
It is currently uncertain how the embryo organizes and directs cellular growth. Several factors are believed to be involved, including the forces at play as cells push and pull on one another [8,41]. Another known factor is that cellular growth can be caused or halted by the introduction of growth factors, chemical signals that are detected by the cells through receptors on the cell surface [51]. Several growth factors exist, signalling cells to react in a particular way, such as undergoing rapid mitosis, halting mitosis altogether, or to migrate to a different area [3]. However, such chemical signals would typically diffuse freely around the cells, lacking the specificity and precision that may be required for proper development. One hypothesis for directed cellular growth is that small extracellular structures called fractones, which are found attached to cells, direct the growth [23]. They do this by capturing specific growth factors with specialized receptors and concentrating them for presentation to the cells the fractone is attached to. In this way, cells with fractones react with specificity to the presence of particular growth factors. Different fractones may target different growth factors, and thus the distribution of specific types of fractones throughout the developing embryos can direct the growth of the entire organism.
Experimental testing of this hypothesis is difficult, however. It is currently impossible to examine how fractones affect development in vivo. For the same reason, it is also impossible to create a time series of development in the same organism. From a biological standpoint, a computational model predicting and matching the observations will support the possibility that fractones are involved in the process of brain morphogenesis.
As biological systems can be extremely complex and since we lack a full understanding of the process of morphogenesis, the model we develop in this paper simplifies the biology, but focuses on the hypothesis that fractones guide development by either promoting or stalling mitosis in specific locations. Several models of cellular growth have been developed in the past [1,53]. Some models work on a tissue-level, examining the growth, motion, and folding of entire sheets or masses of cells [8,10,11,19,52], while others work on a cell-level, examining the dynamics of individual cells and how they interact with each other [36,43], and still others work at a molecular-level and consider the chemical signals that cells use to communicate with each other [51]. As may be expected, some recent models are multi-scale, applying the mechanics of these different scales together [6,7,15,49]. The drawback of modeling on multiple scales is the sheer size and complexity of the system that arises. Tissues can be comprised of millions of cells (the human body has on the order of 10 13 cells [5]), and if each one of those cells is controlled by chemical signals that must be captured by a receptor, then the model becomes computationally intensive. The model presented in this paper is multi-scale, considering diffusion dynamics of chemical growth factors and their capture by fractones. The captured growth factor then triggers behavior changes in the capturing cell, directing a cell to either undergo mitosis immediately or halt mitosis. Replicating cells push their neighbors according to an algorithm, evolving the shape of the cell mass. While extending the model to the tissue-level may seem natural, we are restricted by the complexity of the system to consider only relatively small masses of cells. To our knowledge, as the fractone hypothesis is relatively recent, no hybrid model incorporating fractones has yet been studied and our work is an extension of [9,42].
The model we use is a controlled hybrid automata [26]. Hybrid models have become more important in recent years as technology has increased to the point where near-instantaneous changes can occur in a system due to digital sensors and switches. Hybrid systems have applications to digital technology and automated computer systems. In general, a hybrid system models continuous behavior until certain conditions are met -modeling, for example, feedback from a sensor or a timer -then an instantaneous discrete change occurs which may alter the continuous dynamics. In our model, diffusion of growth factor and capture by fractones is a continuous process, while the growth of cells is a discrete process, occurring instantaneously by assumption, and causing an immediate change to the geometry of the cell mass due to the creation of cells and the subsequent movement and pushing of other nearby cells. The position of the cells blocks diffusion of growth factors, thus each discrete change alters the diffusion mechanics. Our model is complicated due to the size of the state space of the discrete-event system. Each arrangement of cells is a different state, and there are therefore an infinite number of states. The hybrid automata system we adapt for our model formalizes a way in which these dynamics can interact with each other.
In addition to the hybrid and discrete dynamics, we incorporate a control parameter which represents the fractones. Thus the fractones are controlled by the user, and where they exist, they disrupt the diffusion mechanics by acting as sinks. Mitosis that occurs at a neutral rate of once every 24 hours, which is not fractone-driven, introduces a stochastic element as well, as the direction of this neutral growth is chosen randomly. For a brief survey of stochastic elements in hybrid models, see [47]. We note that this random choice of direction is an assumption in the model, as even this non-fractone directed growth is likely controlled in some way in the developing organism.
Our contribution is the development of an innovative model for morphogenesis based on a fractone's governing principle. This demonstrates the applicability of hybrid control systems to complex biological processes.
In this paper, we will first briefly describe the biology behind the hypothesis of fractone-directed growth, and how we hope to model it mathematically. We then formally introduce hybrid automata and formulate our model as a variation of a controlled hybrid automata. We conclude with a discussion of controllability and optimization of the model.

2.
Fractones and morphogenesis. The exact mechanisms that drive morphogenesis are not well known. Indeed, neurogenesis in adult brains was not believed to occur until relatively recently in 1992 [16,22] when specific sites of cellular mitosis in the adult brain, called niches, were found. Recent research has discovered extracellular structures in these niches [33,34], and it has been proposed that these structures control the neurogenesis in the adult brain as well as during early development [23,32]. Named fractones for their branching shape, reminiscent of fractals, these structures are believed to absorb chemical signals called growth factors, and concentrate it for the cells [13]. The location and behavior of fractones therefore directs morphogenesis. See Fig. 1 of the fourth ventricle of an adult rat brain, surrounded by multiple fractones in the subependymal layer.
Fractones were first discovered in 2002 by Dr. F. Mercier, Dr. J. Kitasako, and Dr. G. Hatton. While studying the ventricles of rat brains, they discovered small extravascular bulbs in the subependymal layer. These bulbs tended to be located near groups of replicating cells. Tests revealed that they were in fact basal lamina structures. The basal lamina is a layer of the extracellular matrix upon which stem cells or epithelial cells are located. The extracellular matrix has traditionally been thought of as more of a support structure for the body, however more recent studies have found that it may play a larger role in controlling development [14,45]. The basal lamina is believed to make a significant contribution to morphogenesis and development early in life, but its contributions during adulthood are not clearly known [4,18,50]. The discovery of these basal lamina structures near replicating cells in the subventricular zones under study by Dr. Mercier, et al. was a significant find, as it suggested a mechanism by which neural stem cell growth was regulated.
The recorded basal lamina bulbs measured approximately 1-4 µm in diameter, and were located in a regular distribution close to the walls of the lateral and third ventricles (see Figures 1 and 3), as well as along the developing spinal cord [33,34]. Under transmission electron microscopy, it was seen that these bulbs have a complex branching morphology and are also connected to several cells in the ependyma and subependymal layer (where most cell growth occurs) on one side and blood vessels on the other via thin stem-like structures. As the basal lamina is thought to contribute to morphogenesis, it was hypothesized that these fractones were directly influencing the cells they formed connections with. In 2007, this hypothesis was supported by experiments that found fractones exist in adult humans and have heparan sulfate proteoglycans (HSPG), which capture growth factors (specifically FGF-2) for presentation to cells. The fractones were found to capture and store these growth factors, and that cell proliferation preferentially occurred near fractones rather than blood vessels. Of important note is that not all fractones were found to have HSPG, suggesting that fractones may be specialized in the types of growth factors they will bind [23]. As a result of this study, 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.
Our implementation of the model focuses on the process of neurulation, which is the early development of the spinal column, and leads to the creation of the brain and rest of the central nervous system. The stages of development in neurulation are common between most vertebrates, and is a critically important part of development.
There are many growth factors that are known and their interactions are complex. Different combinations of growth factors can signal different reactions in the cells. Cells typically detect and capture growth factors through various specialized receptors on the exterior of the cell. The specific types of growth factor present will signal cells to replicate, to stop replicating, to migrate, or even to kill themselves (a process called apoptosis). Growth factors may also be responsible for directing stem cells to differentiate into specific types of cells. All of these actions are critical to morphogenesis as it allows cells to create the structures of the body. Growth factor is produced from many sources, including the extracellular matrix and meninges, a sheath-like covering for the brain which has more recently been theorized to play a significant role in development [12]. In our model, the meninges surrounds the cell mass and acts as the primary source of growth factor, producing it at a rate determined by a chosen function. The timing and location of growth factor production offers some determination as to what actions will occur and where, but fractones allow a much finer control than relying on simple diffusion of growth factors.
Hypothesis. It is hypothesized that cellular growth, both in embryos and mature organisms, is guided by the creation and location of fractones. Fractones bind and concentrate specific chemical growth factors. They then present that growth factor to neighboring cells, causing the cells to react according to the type of growth factor presented. In this way, fractones that bind mitosis-promoting growth factors can accelerate growth where they are present, and fractones that bind mitosis-inhibiting growth factors can block growth where they are present.
To simplify the model, only two growth factors are explicitly introduced in our model, with a third implicitly present. The first growth factor represents a positive growth factor, which promotes what we will term accelerated growth every 6 hours. The second growth factor represents a negative growth factor and will inhibit growth. Sufficient quantities of the second growth factor absorbed by a fractone will halt all mitosis for cells adjacent to the fractone. In our model, the two growth factors will diffuse through the system and be gathered by the fractones. The third growth factor is not explicitly modeled, but assumed to be present in the system and absorbed through the receptors directly on the cell, as opposed to absorption by the fractones. The third growth factor is assumed to be responsible for causing cell growth at a slower rate than the positive growth factor. Thus without fractones, cells will still replicate, but will do so at a much slower speed. We will refer to this slower growth as neutral growth in order to reduce the confusion with the accelerated growth directed by fractones.
Studies have since found fractones in many areas of the brain and neural tube, often clustered in areas of high mitotic activity. However, fractones can also be found in areas where cell replication seldom occurs. This lends credence to the theory that growth factors are specialized; some fractones may target growth factors which stop mitosis from occurring. The type of fractone present is therefore just as important as where and when fractones are found during development. In our model, two fractone types will be used since there are two growth factors we explicitly track. The positive fractone will absorb positive growth factor at a much higher rate than negative growth factor, and negative fractones will do just the opposite. It is important to note that in reality, many growth factors are present during development and several types of fractones with different sensitivities may exist. As such our model can be easily extended to incorporate more types of growth factors and fractones.
3.1. Hybrid automata. A hybrid automata has sets of continuous variables and discrete variables. In general, the continuous dynamics evolve the continuous state like a typical continuous dynamical system. This occurs until specific conditions are met. At that point, an instantaneous change occurs in the discrete state. This discrete switch may alter the continuous state, and the continuous dynamics may depend on the discrete state.  [29]). A hybrid automata is defined as a collection H = (Q, X, f, Init, D, E, G, R) where Q is a finite set of discrete variables, and by Q, we denote the set of values these variables can take. X is a finite set of continuous realvalued variables (although any smooth manifold can be used). We denote the set of valuations of n such variables X = R n . An element (q, x) ∈ Q × X is called a state of H. The vector field f : Q × X → T X describes the evolution of the continuous dynamics, and we assume for all q ∈ Q that f (q, ·) is globally Lipschitz continuous. Init ⊆ Q × X is a set of initial discrete and continuous states. D : Q → P (X) is a domain (P (X) denotes the set of all subsets of X), for q ∈ Q, D(q) is the subset of X in which the continuous evolutionẋ = f (q, x) occurs. E ⊆ Q × Q is the set of edges. The edge (q i , q j ) ∈ Q × Q represents the instantaneous change in discrete state from q i to q j . Note that not every pair of discrete states will be in E, as the dynamics of the system may not allow every discrete state to be reached by every other discrete state. G : E → P (X) represents the guard conditions for each edge, the subset of X which will cause a switch in the discrete state, along the given edge. Finally, R : E × X → P (X) is the reset map for each edge. When a discrete switch occurs along edge E, X may change, causing a discontinuous jump in the continuous dynamics.
Examples of hybrid automata can be found in, [20,27,29]. In the description of the hybrid automata, H, there are very few restrictions on the discrete mechanics described by the guard conditions, edges, and domains. It is generally incumbent on the model creator to ensure that the hybrid automata will operate as intended. The following two definitions help to account for a properly working hybrid automata.
Definition 3.2 (See [29]). A hybrid time trajectory, τ , is a collection of intervals for continuous growth: τ = {I i } N i=0 such that: The idea is that a hybrid time trajectory gives us the time period in which we are in a sequence of discrete states. Note that discrete changes are assumed to occur at τ i = τ i+1 . We further define τ = {0, 1, . . . , N }, the indices of the intervals in τ . And we denote |τ | = i∈ τ (τ i − τ i ). [29]). An execution of a hybrid automaton, H, is a collection χ = (τ,q,x), where τ is a hybrid time trajectory,q : τ → Q, andx = {x i (·) : i ∈ τ } is a collection of differentiable maps x i :

Definition 3.3 (See
Thus an execution, if one exists, describes a proper path in the system where we start in a discrete state, the continuous dynamics flow according to f until we hit a guard condition at the boundary of the domain of the discrete state, and this immediately changes our discrete state and our continuous state changes (or in some cases remains the same) so that we are in the domain of the new discrete state, and the continuous dynamic resumes (although potentially affected by the discrete state) until we hit another guard condition. Determining whether a hybrid automata has executions for any given initial conditions, and whether those executions are unique, is one important area of study [27]. Continuity and stability of hybrid systems can also be analyzed [28]. While typically hybrid automata are sensitive to changes in initial condition, especially as they become more complicated in modeling real-world systems, it is possible for hybrid automata to have some sense of continuity. Briefly stated, a hybrid system is continuous if two executions which are close remain within a given distance of each other, where this distance is determined by a given metric for Q × X and the difference between the total time |τ | of the executions. Definition 3.4 (See [29]). An execution is finite if τ is a finite sequence ending with a compact interval. We denote by ε * H (q 0 , x 0 ) the set of all finite executions. An execution is infinite if either τ is an infinite sequence, or |τ | = ∞. We denote by ε ∞ H (q 0 , x 0 ) the set of all infinite executions. An execution is zeno if τ is infinite and |τ | < ∞. We denote ε H (q 0 , x 0 ) as the set of all execution of H with initial condition (q 0 , x 0 ).
Finally, we introduce the notion of reachable sets.
Definition 3.5 (See [29]). We define the set of all reachable states by H as x)} Analyzing the set of reachable states is important for the next question of controllability. Several papers cited above deal briefly with the question of controllability in hybrid automata. There is also similar work done on different types of hybrid systems [54], many of which can be rewritten as hybrid automata [26], and some work on controllability for general hybrid systems [55] also exists. The general algorithms given, however, are too computationally intensive for a system as complicated as our model and we will need to develop ad-hoc techniques.

3.2.
Fractone hybrid automata. In this section, based on our assumptions from section 2 we model the role of the fractones in the morphogenesis process as a hybrid automata. Morphogenesis is a time-invariant process, therefore we can assume without loss of generality that our hybrid time trajectories satisfy τ 0 = 0.
3.2.1. Discrete-event system. The discrete-event system of our model describes the asynchronous events triggered by mitosis. The state space corresponding to the discrete-event system represents the physical and chronological distribution of cells, meningeal cells, and fractones in the ambient space. Although cells can take a variety of shapes and sizes, and to a degree tend to be malleable, we simplify our model by assuming cells are ellipsoids of the same volume, and do not deform after creation. Our definition of cells as ellipsoids of fixed volume but variable semi-axis length has been used in previous models [24,37,46,56].
Definition 3.6. We define a cell body as an ellipsoid of fixed volume V and semi-axis lengths r 1 , r 2 , r 3 , with s ≤ r i ≤ S, i ∈ {1, 2, 3}, for some constants V, s, S ∈ R. For our model, we use an average cell size, relying on measurements taken of developing rat embryos showing the neuron soma size to be approximately 72 − 81 µm 2 in size [31] and averaged 9 µm cell diameters of yeast cells [56], and set s = 3 µm, S = 6 µm, and V = 4 3 π(4.5) 3 µm 3 . Although mitosis in our model occurs instantaneously, we must account for the time it takes a mother cell to split into two daughter cells before the new cells can undergo mitosis. This time for cell division in an embryo varies depending on the type of cell, age, and species. In a developing rat embryo, this time can vary between 3 to 14 hours [2]. For our model, we use an approximate 6 hour cell cycle as an assumption. To each cell body, we introduce a parameter called cell-birth to account for the age of the cell and determine when a new cell is ready to replicate. It is measured in minutes. It is important to note that when mitosis occurs, the mother cell splits into two daughter cells, thus the parent cell ceases to exist and two daughter cells are created, both with cell-birth parameter set to the time of cytokinesis (when a single cell is divided to form two daughter cells).
Definition 3.7. To a given cell body c, we associate a cell-birth parameter λ c to track the age of the cell. A cell is defined as a pair {c, λ c } where λ c is the birth time of the cell with respect to τ 0 = 0, measured in minutes. A cell-birth parameter can be negative to account for possible age differences of the cells at the beginning of the process.
In the sequel, we restrict the possible sets of cells to those that make physical sense; cells should not overlap with each other in space and a mass of cell should be connected. To determine how far apart cells may be in an acceptable set of cells we introduce a new parameter . Distances greater than 2 between cell boundaries will cause holes to appear which we would like to avoid. In our model, we use = 4.5 µm, the same length as the radius of a spherical cell.
,··· ,n be a set of n cells and B = {c i } i=1,··· ,n the corresponding set of cell bodies with B ⊂ R 3 . To each cell body c i , with semiaxis lengths r 1i , r 2i , r 3i we assign a concentric ellipsoid,ĉ i , with semi-axis lengths r 1i + , r 2i + , and r 3i + respectively. LetB = iĉ i . IfB is a compact connected space of R 3 and c i ∩ c j = ∅, ∀i, j ∈ {1, · · · , n}, i = j, then C is said to be admissible. Figure  4. Notice that it is a purely geometrical notion with the age of a cell being irrelevant. Meninges will play an important role in our model, and for clarity we will always refer to meningeal cells as such and embryonic cells, as described previously, as just cells. Meningeal cells are located on the exterior of a developing set of cells, and are the source of growth factors in our model. To simplify the dynamics of the model, meningeal cells undergo no mitosis and their age is therefore unimportant. Definition 3.9. We assume that meningeal cells are of uniform volume and spherical shape. More precisely, the meningeal cell c m centered at (x 0 , y 0 , z 0 ) ∈ R 3 is defined as the closed ball of radius 2 centered at (x 0 , y 0 , z 0 ). A meningeal cell in our model has a radius of 2.25 µm, with a volume one-eighth that of a spherical cell.

Examples of admissible and inadmissible sets of cell bodies can be found in
As with sets of cell bodies, we restrict the location of meningeal cells to an admissible set.
Definition 3.10. Given an admissible set of cells, C, and a set of meningeal cells, C m = {c mp }, which we will call the meninges, C m is admissible if the center of each meningeal cell lies on the boundary, ∂B, ofB, and c mi ∩ c mj = ∅, ∀i, j, i = j. By definition, a meningeal cell cannot intersect a cell.
The developing organism is primarily a group of cells surrounded by meninges, we therefore define their combination. Finally, we introduce the fractones to our discrete-event system. Fractones have a tree-like shape when viewed with electron microscopy, however, for this paper, we do not study the fractone's specific geometry, its receptors, or its means of attachment to cells, but only its specific hypothesized role in development. We thus model fractones as small spherical structures found next to cells. Fractones are assumed to be closed balls with radius one-ninth that of a spherical cell, i.e. they have a diameter of 1 µm, which is on the smaller end of the 1-4 µm diameter found experimentally. It also guarantees that they will never intersect with meningeal cells.
where r = 1 9 3 3 4π V = 0.5. Given an admissible set of cells, C, a set of fractones, F = {f i } i=1,··· ,p , is admissible if every fractone is tangent to at least one cell body and f i ∩ f j = ∅, ∀i, j ∈ {1, · · · , p}, i = j. A fractone f is tangent to a cell body c if f ∩ c is formed of exactly one point.
The combination of cells, meningeal cells, and fractones are combined to form the state space of our discrete-event system, which we call a biological structure. As mentioned in the introduction, we only consider two types of fractones to reflect accelerated and prevented growth. Definition 3.13. We define a biological structure q as a triple, {M, F + , F − }, as a cell mass M , paired with two admissible sets of fractones, F + and F − , which represent the positive and negative fractones, respectively, in the system. Q, the state space of the discrete dynamic of our model, is the set of all possible biological structures in the ambient space.
This differs from traditional hybrid automata studied in the literature since the state space, Q, is not finite. A cell can be located anywhere in the ambient space, so there are uncountably infinite positions for a cell to exist. The same is true for meningeal cells and fractones, even though all elements of the biological structure have restrictions on placement in order to be admissible. Notice that the cell-birth parameter also introduce an infinite dimension.
In the literature of hybrid automata, the topology of the discrete state space is taken as the discrete topology generated by the metric d D (q, q ) = 1 if q = q and d D (q, q ) = 0 if q = q for q, q ∈ Q. However, since the state space corresponding to the discrete-event system is not finite nor countable in our case, and due to the specificity of our application we introduce a different notion of distance. The motivation comes from the fact that in comparing two different organisms of the same species and age, we expect some differences between them -the exact length and width of limbs, for example -but the general shapes ought to be similar; limbs should not be growing from different places or missing altogether. Our metric is based off of the Hausdorff distance [44] and the Minkowski metric as defined on the Minkowski spacetime [35]. Indeed, since the definition of biological structures includes cell-birth parameters, we need a notion of spacetime interval.
Definition 3.14. Consider two biological structures, ,··· ,na , a = 1, 2 the corresponding admissible sets of cells. We define the geometric distance between q 1 , q 2 by d g (q 1 , , where d H denotes the Hausdorff distance for the Euclidean metric and B a = {c i } i=1,··· ,na , a = 1, 2. To account for a possible age difference of the cells between two biological structures, we also introduce the age distance d a . Notice that this temporal distance also involves a spatial component to account for a measure of the age difference with respect to all cells. We denote the Euclidean distance between the centers of two cell bodies c 1 For our model, we use κ = 4.5, which is chosen so that a 1 minute difference in age is equivalent to spherical cells being a radius length off-set. Its unit is microns/minutes, and it represents an arbitrary standard velocity chosen to express how far apart two cells are. Finally, we define the distance between q 1 and q 2 as d(q 1 , q 1 ) = d g (q 1 , q 2 )+d a (q 1 , q 2 ). See Figures 5 and 6 for examples.  Proof. It follows directly from the fact that our distance is the addition of d g which is a metric by definition and of d a which is defined as the Hausdorff distance of a metric defined by the Euclidean metric between cell body's centers and a multiple of the absolute value of the corresponding ages. It is a pseudo-metric because two distinct q 1 , q 2 ∈ Q with the same fractones and cells but different meninges imply d(q 1 , q 2 ) = 0 even though they are not equal.
Note that the pseudo-distance would be a distance if we introduce an equivalence relation to identify cell masses that differs only by their meninges.
3.2.2. Continuous-event system. Growth factors diffuse around cells and meningeal cells and cannot pass through them, while being absorbed by fractones. We define the diffusion mechanics for two growth factors, but the model can be extended naturally to more than two growth factors.
For a given q ∈ Q, growth factors will diffuse only through a space which we will call the diffusion space. This diffusion space extends a certain radius around the exterior of the set of cells or more depending on the region of interest, as growth factors could realistically be assumed capable of diffusing across cavities in a developing embryo. Cells and meningeal cells are removed from the diffusion space as we assume there is no diffusion of our two growth factors across the cell membranes. The particular geometry of the diffusion space extends out so that the center of the meningeal cells lie on its boundary. This is done because the meninges acts as an external sheathe for the developing organism, and growth factor would not be expected to diffuse across it. Examples of biological structures and diffusion spaces can be found in Figure 7. Note that the diffusion space changes as the biological structure evolves.
Definition 3.15. The diffusion space corresponding to a given biological structure q, denoted by D q is the space in which the growth factors diffuse for this biological structure. It is determine by q and the region under study. In practice, the initial diffusion space will be determined by the experimental maps. For traditional hybrid automata, the continuous dynamic is defined on a finite set of continuous real valued variables. In our case, we need to expand this notion since our continuous dynamic must be defined on an infinite dimensional space. Moreover, that infinite dimensional space varies with the discrete dynamic which creates an even more intricate model. For a fixed biological structure q, we introduce the set of functions X q : D q ⊂ R 3 → R ≥0 × R ≥0 that we denote by X q . We have for x ∈ D q that X q (x) = (X 1,q (x), X 2,q (x)) where X 1,q (x) (respectively X 2,q (x)) represents the concentrations of positive (respectively negative) growth factors at x = (x 1 , x 2 , x 3 ) ∈ R 3 . The units of X q (x) are arbitrary molarity units, which we will term concentration units. A state of our hybrid automata H is given by a pair (q, X q ) ∈ Q × X q .
The evolution of the continuous dynamic is determined by the diffusion of the growth factor and their production. In the sequel, for a given q ∈ Q we will denote the concentration of growth factors at time t and position x in the diffusion space D q by X q (t)(x) = X q (t, x). The continuous dynamics is given by the following equation: where∆ represents a perturbed diffusion with (κ 1 , κ 2 ) the diffusion constants of the positive and negative growth factors, and g i (t, x, X i,q ) accounts for the production of positive and negative growth factor. Note that in the absence of fractones and meningeal cells (which imply no production of growth factor), diffusion occurs in the diffusion space according to the heat equation. When fractones are present, diffusion is perturbed. The fractones act as sinks, with growth factor only flowing into fractones, and never out. This means that the spatial distribution of fractones controls the diffusion process and by extension the evolution of the biological structure. To describe the effect of the fractones on growth, we introduce a control function u : D q → {0, 1} × {0, 1} that we assume measurable bounded. The role of the control function is to switch on and off the existence of fractones in the diffusion space. More precisely, let q = (M, F + , F − ) be a biological structure and f + be a positive fractone centered at (x 0 , y 0 , z 0 ), i.e. f + = B r (x 0 , y 0 , z 0 ) ∈ F + . For every x ∈ B r (x 0 , y 0 , z 0 ) the control function takes the value u(x) = (1, 0) to indicate the existence of the positive fractone. More generally, we have that u( The perturbed diffusion then is modeled accordingly to the control function u to capture the role of the fractones as sinks in the diffusion process. Where fractones exist, the diffusion space appears empty of growth factor when calculating the diffusion mechanics. Fractones do store growth factors -they merely make them unavailable for diffusion. Thus the control u has an influence on the diffusion mechanics. The perturbed diffusion can be made explicit using a discretization of the ambient space, equation (1) can then be written as an affine control system.
The production of growth factor is represented as a reaction term in our equation with the restriction that g i (t, x, X i,q ) ≥ 0. To ensure existence of a solution for the continuous dynamics we assume g = (g 1 , g 2 ) is continuously differentiable with respect to all of its variables. Growth factor is introduced into the diffusion space through the meninges, it is indeed assumed to be produced by the meninges and is added to the diffusion space near the border of any meningeal cells. For our simulations, the volume in which the growth factor is added at the time of production corresponding to a biological structure q, is defined by Disk m = B /2+ g (x 0 , y 0 , z 0 ) ∩ D q for each meningeal cell, c m , centered at (x 0 , y 0 , z 0 ) where g = 0.5 (an arbitrary choice). At this stage there are only speculations about how the growth factor production might vary and no conclusive experimental results exist. It is expected that the mathematical model can be used to try many different possibilities for growth factor's production to guide the experimentations and improve our understanding of this phenomenon. In our model, positive and negative growth factor production may be different and they can depend on time, position and on the actual concentrations (for instance if there is a high concentration of a given growth factor at the meninges the system might prevent its production). For the simulations presented in section 5, we assume an equal production of growth factor by each meninges cell. We have conducted some preliminaries testing with constant functions, sinusoidal functions, and square waves, which are all assumptive of how growth factors might be produced by the meninges.
To initialize the model, we need a starting point that includes the initial biological structure and the initial distribution of growth factors. In practice, we obtain these initial conditions from biological experimentation.
Definition 3.17. The initial conditions are dependent on the age of the embryo and location under study. The initial conditions for the discrete dynamics and continuous dynamics are given in Init ⊆ {(q, X q ); q ∈ Q}.
As an example let us consider an initial condition starting at time t = 0 with one cell centered at x = (x 1 , x 2 , x 3 ) = (0, 0, 0), no meninges, a distribution of the first growth factor concentration which grows from 0 to 20 concentration units as we move from the boundary of the cell to the boundary of the diffusion space, and a distribution of the second growth factor which is 0 everywhere in the diffusion space, see Figure 8. The initial conditions are then given by Init = {q, (X 1,q , X 2,q )}, where  The interplay between the discrete and continuous dynamics is modeled through several mathematical structures: the domains, edges, guard conditions, and reset maps.
Definition 3.18. The domain, Dom, assigns to each q ∈ Q all of the growth factor distributions that could occur without a mitosis event. We have Dom : Q → X × R ≥0 , where X denotes the set of all X q . For a q ∈ Q, Dom(q) maps to the set of all growth factor distributions of our two growth factors and all times in which cellular growth has not occurred.
Each mitosis event changes the biological structure, and the possible changes that could physically occur are represented in the set of edges, E ⊆ Q × Q. Definition 3.19. An element (q 1 , q 2 ) ∈ E, the set of edges, represents the instantaneous change in the discrete dynamic from biological structure q 1 to biological structure q 2 . Note in our system, not every element of Q is connected to every other element. The existence of an edge between two biological structures determines the rules for growth, and is represented in our model by our control and the algorithm for pushing cells.
To determine which edge is traversed and when, each edge is assigned a guard condition. For an edge (q 1 , q 2 ), when this guard condition is met, growth occurs instantaneously and we instantly switch from biological structure q 1 to the new biological structure, q 2 .
Definition 3.20. There are two sorts of guard conditions, one for the accelerated growth and one for the neutral growth. The guard conditions for accelerated growth represent those growth factor distributions where a fractone has absorbed sufficient positive growth factor, and the guard conditions for neutral growth represent when sufficient time has passed for a cell to undergo neutral growth. Notice that the negative fractones inhibit growth and therefore prevent a guard conditions to be met. The guard conditions are given by G : E → X × X × R ≥0 . For a given edge (q 1 , q 2 ), G ((q 1 , q 2 )) maps to the set of all growth factor distributions of our two growth factors and times at which mitosis, either accelerated or neutral, would be triggered, and lead to the biological structure deforming from q 1 to q 2 . See Figure 9 for a visual representation of how the domains, edges, and guard conditions interact. Figure 9. The interaction of the domain, guard conditions, and edges of our model. Left to right, we begin in the domain of biological structure q 1 . The guard condition for edge (q 1 , q 2 ) is then met when the positive fractone captures 100 concentration units of growth factor. This triggers growth and we instantaneously switch discrete states to biological structure q 2 .
When growth occurs and new cells are created, growth factor will be redistributed in the system due to the physical forces related to the motion and growth of the new cells. This is modeled in our system by a deformation algorithm, which is represented mathematically by the reset map, see below. The guard condition displayed in Figure 9 is typical, and when it is met, mitosis occurs, causing an instantaneous change in the state of the discrete-event system (along some edge). The edges define the growth algorithm. As has been stated, we do not model the actual growth of the daughter cells -all daughter cells are formed at full adult size immediately. When a cell undergoes mitosis, the daughter cells are formed along the axis of the tangent vector to the surface of the ellipsoid at the point where the fractone is tangent to the cell. In this way, the fractones determine the general direction of growth. The cell bodies of daughter cells are separated by a fixed distance d m along this tangent axis when placed. The value for this distance is motivated by the fact that we do not want to create holes greater than the size of a single cell, which would be d m < 2r where r is the radius of a cell. In practice, we would use something close to d m = r (for our numerical simulations, we use r = 4.5 and d m = 3). The new cells have the same geometry and orientation as the parent cell by simplifying assumption, and when it is created it may force pushing of neighboring cells that may cascade into a deformation of the entire cell set. The directions of each push is determined by the direction of the tangent vector to the surface of the ellipsoid being pushed at the point of collision. A simple example is shown in Figure 10. The growth that occurs pushes growth factor out of the way according to the reset map. Let V gf be the volume formed by the union of the new daughter cells and the elliptic cylinder, V cyl , with faces formed by taking the planar cross-sections of the daughter cells, perpendicular to the aforementioned tangent vector, through their centers. You can think of V gf then as the volume formed by starting with both daughter cells centered at the location of the mother cell and then sweeping out to their proper positions along the direction of the tangent vector. All of the growth factor in V gf is evenly redistributed in a volume along the surface of the daughter cells, depending on how the growth occurs. For example, in Figure 10, only one of the daughter cell moves from the position of the mother cell, so all of the redistributed growth factor is distributed around that displaced cell, and none around the daughter cell that remained in same location as the mother cell. In particular for this example, we set a parameter gf , and define V out , the set of all points on the exterior of the daughter cell a distance gf from its surface. Then the volume in which the growth factor is redistributed is V out \V cyl . A similar scheme is used for growths in which both daughter cells are displaced from the position of the mother cell, one simply finds an analogous V out for both cells.
Definition 3.21. For each edge, the reset map determines how the growth factor distribution changes when the guard condition is met. Thus, although we call the diffusion the continuous dynamic, each mitosis event may cause a discontinuous jump in the growth factor distribution over time. The reset map is given by R : E × X → X . As mentioned, upon triggering a growth along edge E, the values of X are changed, representing the pushing of growth factor due to the growth. For mitosis occurring at time t m along edge (q 1 , q 2 ), X q2 (t m , x) = R ((q 1 , q 2 ), X q1 (t − m , x)), where t − m is the limit as t approaches t m from the left.
It is important that the initial conditions are within the domain of the initial biological structure, and that our reset maps similarly always redistribute the growth factors so that we are within the domain of the new biological structure. This way, the system is always within Dom(q 1 ) for some q 1 ∈ Q, and when it meets a guard condition for some edge (q 1 , q 2 ) ∈ E, which triggers mitosis and a discrete switch to the new biological structure, q 2 , the reset map ensures the system is now in Dom(q 2 ). Let us consider an example represented in Figure 11. We have a biological structure, q with two growth factors, cell bodies c 1 , c 2 , c 3 , cell-birth parameters λ c1 , λ c2 , λ c3 respectively, positive fractones f + 1 , f + 2 , negative fractone f − 3 , and diffusion space D q . Growth will occur when one of the positive fractones absorbs enough positive growth factor or the time reaches the neutral growth time of one of the cells. Cell c 3 will only undergo neutral growth if the negative growth factor absorbed by f − 3 does not exceed the threshold. The thresholds are arbitrarily chosen, as our growth factor concentrations are in arbitrary units. These thresholds represent the level of growth factor concentrations at which the cells react to the growth factor, either undergoing mitosis for positive growth factor or halting all mitosis for negative growth factor. For an accelerated growth the threshold is taken as τ 1 = 100, and the negative growth threshold as τ 2 = 80. Finally, the neutral growth time delay is fixed to 360 minutes. For this biological structure, we have that where V represented the volume of the fractone.
A mitosis occurring for cell c 1 would be represented by an edge, e 1 , with guard condition: 4. Controllability and optimization. We make the following definitions to describe the evolution of the system over time.
Notice that given a hybrid time trajectory τ , initial conditions (q(0), X 0 (0, x)), and a collection of admissible controls as in definition 4.1 an execution for our system always exists as long as we permit the diffusion space to grow with the evolution of the biological structure. However, because of the random nature of the neutral growth, we have that an execution is not unique. Definition 4.2. Let χ = (τ,q, X, u) be a finite execution and a fixed initial condition (q 0 , X 0 ). We define its end-point set Λ H (q 0 , X 0 , χ) as the set of all possible end-points (q(N ), X(τ N )). We have that Λ H (q 0 , X 0 , χ) ⊆ Q × X .
Since the random neutral growth tends to complicate the system by adding a stochastic element and many end-points for a given control, we will often consider a simpler system by removing the random growth mechanic. In this simplified system the growth of cells occurs only when a fractone absorbs sufficient growth factor, and the evolution of the system is now completely determined by the initial conditions and control used (we assume the production of growth factor has been fixed). This means that Λ H (q 0 , X 0 , χ) is in this case a single end-point.
The set of all biological structures that can be created from a fixed initial biological structure during a morphogenic process taking place over a fixed duration is defined as the evolution set. The evolution set of H from q 0 is given by Evol q0 H = q0 T ≥0 Evol H (T ). A zeno execution, named after the philosopher Zeno and his well-known paradoxes, is one that takes finite time but has an infinite number of switches. Authors in [20,21] have examined zeno executions, specifically examining ways to extend the zeno executions beyond their zeno time (the finite time of the execution). This is important because in simulating hybrid automata, zeno executions create problems when they begin to switch very quickly -lots of computations occur due to the switching, but no real progress is made since time barely moves. Methods for breaking the zeno behavior have been proposed in order to extend a zeno execution to an infinite execution. Given the nature of our problem and our hypothesis, it is natural to suspect that mitosis occurs in a discrete way. Indeed, it can be proved that the fractone hybrid control system H does not have zeno executions under a biological relevant condition on the production of growth factors.
Proposition 2. Zeno executions do not exist for the fractone hybrid control system provided that the production of growth factor is bounded over hybrid time trajectories τ with N < ∞.
Proof. Since cells, meningeal cells, and fractones must exist in an ambient space, which is assumed to be of finite volume, there can only be a finite number of each of them at any time. Assume χ is an execution with N < ∞. Since we have a finite number of meningeal cells, and using the hypothesis that the production of growth factor is bounded we have that the maximum amount of growth factor that could ever be in the system is finite. Let this maximum amount be G. To undergo accelerated mitosis, a fractone must absorb at least τ 1 positive growth factor, and this absorbed growth factor is removed from the system afterwards. Thus at maximum, there are G/τ 1 accelerated mitosis events. In addition, for a growth time delay, t d > 0, there are at maximum, τ /t d neutral growths. Both of these numbers must be finite since G and τ are finite, thus there are only a finite number of discrete events that can occur in finite time. Therefore Zeno executions are not possible.

4.1.
Controllability. Our on-going research revolves around this question of controllability as we hope to simulate the development of embryos over time using the hypothesis that fractones guide growth. Given experimentally-determined maps of developing embryos at different times, we hope to be able to find the control that takes us from one map to another. For example, if we know the location of cells in a portion of a developing rat embryo seven days after fertilization and eight days after fertilization, we hope to be able to find the distribution of fractones (the control) such that our simulation takes us from day seven's cell map to a reasonable approximation of day eight's cell map. Finding such a control, and verifying its biological plausibility, would lend credence to the hypothesis that the development is guided by the fractones.
Definition 4.4. Let q be a fixed biological structure. We say that the system is controllable to q from q 0 in time T if there exists X ∈ X q such that (q, X) ∈ Evol q0 H (T ). Notice that since we are considering biological growth, we must grow at a particular pace and reach our target at a particular time, with a certain degree of accuracy. This accuracy, however, need not be perfect. Problem 1. Let δ > 0. Given an initial biological structure q 0 ; a target biological structure, q 1 ; an initial growth factor distribution, X 0 ∈ X q0 ; a fixed time T > 0; can we find an execution with |τ | = T , and such that d(q(N ), q 1 ) ≤ δ?
Problem 1 is a very complex and restrictive problem since we need to take into account the fact that we impose the initial growth factor distribution (it is therefore a constrained controllability problem) and require the final configuration of meningeal cells, and birth-parameters to be close to the desired biological structure. However, it leaves completely free the final growth factors concentrations. Moreover, the random nature of the neutral growth makes it very difficult to analyze. We therefore state a sub-problem by relaxing some of the constraints.
Problem 2. Let δ > 0 and let us remove the neutral growth. Given an initial biological structure q 0 ; a target biological structure, q 1 ; an initial growth factor distribution, X 0 ∈ X q0 ; a fixed time T > 0; can we find an execution with |τ | = T , and such that d g (q(N ), q 1 ) ≤ δ?
The main difference between problem 1 and problem 2 is that we impose a control over the final set of cell bodies only rather than on the entire biological structure, and of course the removal of the stochastic aspect of our system. The literature on hybrid systems has limited results on controllability of a system as large and complicated as ours. Most papers focus on classes of hybrid control systems that are not applicable to ours [21,20,54] or look at general hybrid control models but propose algorithms to determine controllability that, while compelling, are in practice infeasible for the complexity of our model [55]. However, for problem 2, given the flexibility on the spatial distribution of fractones we can prove the following result.
Proposition 3. Given an initial biological structure q 0 ; a target biological structure, q f with the corresponding sets of cell bodies B i satisfying B 0 ⊆ B f ; an initial growth factor distribution, X 0 ∈ X q0 , then there exists a finite execution χ = (τ, q, X, u) such that X(0) = X 0 with d g (q(N ), q f ) ≤ 13.5 provided that neutral growth is removed from the system dynamics.
From a practical point of view, we will allow the cells to be a little more than one cell-width off from the target. Since a spherical cell in our model is 9 µm in diameter, the δ used in the controllability proof above is biologically relevant.
Proof. To allow for a maximal control over the evolution of the biological structure, the idea is to construct an execution with the existence of only one positive fractone at a time. The first step in the construction is to determine the number of cell bodies to be created and their locations. This is done by adding cell bodies to B 0 spaced a distance d m = 4.5 µm apart from a neighboring cell until the new mass of cell bodies, denotedB, contains B f and such that d g (B, B f ) ≤ δ. It is straightforward to prove that there exists an execution that will produce the mass of cell bodiesB. Indeed, the biological structure q 1 is determined by choosing a cell body in B neighboring B 0 and place a positive fractone on the exterior of B 0 such that the triggered growth will occur in that direction, producing a cell body in the location we chose. We can always chose a production of growth factor g 1 to be such that the continuous dynamics guarantees that eventually the positive fractone absorbs enough positive growth factor to trigger accelerated growth. After growth has occurred, remove the fractone and repeat the previous step by creating a new biological structure by placing a positive fractone in an appropriate location to trigger growth as desired. Repeat until the final target is reached.
This new mass of cell can always be grown simply using the above single-fractone algorithm since every the tiling was created progressively from the initial B 0 . We note that different sets of cell bodies can be created, all within geometric distance 13.5 µm. Thus solutions to the controllability problem exist, and they are not necessarily unique. Growth in this single-fractone fashion will never be blocked and never cause pushing because of our choice of fractone position and our modification of the target structure. Modifying a target structure that contains ellipsoids of different shapes is a straightforward adaptation of this process. Proposition 3 provides only a partial answer to problem 2. Indeed, this algorithm provides a solution that is very slow, as it only allows one fractone at a time, and we must wait for growth to occur. However, since the fractone is always placed facing the exterior of the cell mass, it is also very close to the meninges, so growth factor reaches it very quickly to trigger the accelerated growth.
We remark that removing the neutral growth from the system is a significant change to the dynamics of the system. However understanding the behavior of the system without neutral growth is a first step to understanding the dynamics with it. The next step would be to evaluate how significantly the results from the simplified system are perturbed by the reintroduction of the neutral growth, and whether the question of controllability we have stated can be altered to compensate, perhaps through relaxing the bound on the distance between the final set of cells and the target set by an acceptable amount.

4.2.
Optimality. Although it is beyond the scope of this paper, the related question of optimality is mathematically interesting and potentially biologically relevant. Biological mechanisms often become quite efficient by some measure through the natural course of evolution (although exceptions abound). Comparison between results for optimal controls of various cost functions defined on the system may reveal general strategies adopted by the developing organism for efficient growth.
Assuming a biological structure is reachable from a given biological structure, we would like to demonstrate the existence of an efficient control based on some biological relevance. This cost function could measure several factors, including time, accuracy, or the number of fractones we use in the system, on the assumption that the creation or destruction of fractones would require energy from the developing organism.
To determine a cost function, we intend to use a computational model to compare with experimental maps, determining controls that achieve the growth that is found in nature. By then analyzing these controls and the distribution of the fractones, we hope to find any patterns or behaviors that suggest a strategy, and build a cost function based on the observations. We implemented a numerical version of our model in MATLAB 2010b. The numerical model uses some simplifications to the original model. For example, structures in the simulation are aligned to a three-dimensional grid, and all cells are the same size and shape -in particular, they are 9µm × 9µm × 9µm cubes. Representing the cells as cubes instead of spheres simplifies the programming of the simulation. The time is discretized into time steps. In the simulations displayed in this paper, we use 200 time steps for each day (thus 1 time step is 7.2 minutes). As in the model, growth factors are created and diffuse through the diffusion space using standard numerical techniques. When sufficient amounts are captured by fractones, they trigger growth according to the growth algorithm. This growth is immediate and instantaneous and alters the diffusion space based on our deformation algorithm. Growth factor diffusion and absorption then continues. We present two simulations here. Both simulations test the flexibility and performance of the model. We assume for both a production of growth factor corresponding to an added units of concentration of 15 positive and 15 negative in the diffusion space adjacent to a meningeal cell every 10 time steps. While this is not a differentiable function, it is easily handled numerically.

5.1.
Folding of the neural tube. The first simulation is inspired by a portion of neurulation, the folding of the neural tube during development. At the start of neurulation, we have a neural plate which folds up to form a U-shape, and then continues to fold to create a hollow cylinder. We simulate the growth from the Ushape to the tube (see Fig. 13). This simulation just demonstrates how the model behaves, and is not a realistic simulation in terms of size or biological mechanisms of actual neurulation. Cells in our model, for example, do not migrate, which is an important part of the neurulation process. Figure 13. Top view of the neurulation simulation. In black, the cells. In grey, the fractones. We start with a U-shape configuration and wish to end with a closed ring of cells.
Closing the U-shape in the model for our choice of parameters takes 1150 time steps (5.75 days). The progression of the growth is shown in Fig. 14. Note that we have removed the slow neutral growth from this simulation. The reason is that since the neutral growth is random, it makes matching our desired shape more complicated. As a consequence, growth occurs only by capture of growth factors by the positive fractones. No negative fractones were used due to the lack of neutral growth, however when neutral growth is reintroduced in future simulations, we intend to examine their impact on controlling the random growth. Placement of the fractones is determined by the choice of the control function, and fractones are redistributed every 50 to 100 time steps in order to guide the creation of the tube. At various times, up to 90 positive fractones and as few as 10 positive fractones are used. We start with 1080 cells and end with 2110 cells, thus there are 1030 mitosis events. Since growth is guided only by the fractones, and growth factor was plentiful, growth occurred approximately every 50 time steps (according to the 6 hour period it takes immature cells to be ready to reproduce). Note the strategy we used to close the tube, first placing fractones to elongate the sides of the U-shape, then closing the tube, and finally rounding out the top of the tube to be symmetrical with the bottom. The execution shown is simple and clearly not unique -it is not difficult to see some refinement could be done to reduce the amount of time necessary to form the tube. While negative growth factor was produced in the system, the absorption rate of the negative growth factor by the positive fractones was relatively small compared the absorption rate of the positive fractones (0.2 compared to 0.9 respectively), so the effect of the negative growth factor was limited.

5.2.
Controllability for a complex shape. In our second example, we offer an example of a initial set of cells, E 0 , and final set of cells, E f , with some complexity and a control that guides the hybrid automata developed in this paper from one to the other. This control solution is once again not unique, but we do have some goals in creating it. First, of course, we want accuracy, so we impose the geometric distance to be less than 12, i.e. d g (E 0 , E f ) ≤ 12 (we will ignore the age distance for now). Second, we want to be biologically-feasible. This means we will not constantly create and destroy fractones; we will try to maintain the same fractones in the system for as long as possible. It is known that fractones do not migrate, and suspected that they have a life duration between 3 to 12 hours. Since this simulation takes place over ∼ 5 days time we expect only a limited number of changes in fractone location. Our target and initial cell configurations are shown in Fig. 15.
We generate the approximation to the final state with a control, u(t, x). At the start of the simulation, u(0, x) describes the initial fractone positions, and only at specific times (shown in Fig. 16) do we set u(t, x) to not follow the usual growth and displacement algorithm for relocating fractones and instead create and remove fractones in specific locations to generate the approximation to the final state.
To create the control that generated this final target configuration of a rabbit, positive fractone positions were set by the user and growth allowed to proceed until the growing cell configuration exceeded the target cell configuration's dimensions. As with the previous simulation, the random neutral growth was removed from the model. The state of the system was saved at every 5 time steps, allowing the user to then look through the evolution of the growth and choose a point in time to alter the fractone distribution (these are the exact times shown in Fig. 16) in order to Figure 14. Progression of the growth of the tube. Shown is the mass of cells at selected times. Cells (black), fractones (light grey), and immature cells (dark grey). Meninges is not shown in the image for clarity -it is present surrounding the mass, however. From t = 0 to t = 210, we simply lengthen the sides of the Ushape. Until t = 560, we continue to extend upwards, but with staggered growth to create the curved edge. Around t = 1000 we have closed the ring, and at t = 1150 we have completed the top of the ring.
avoid exceeding the target cell configuration's bounds, and produce a reasonable approximation of the rabbit. Note that early on, a large amount of time (200 time steps, or 1 day) passes before fractones are altered, but as more detailed growth is required, the time between changes shrinks to about 50 time steps (6 hours). The final cell configuration and the target cell configuration have a geometric distance of 12. The oldest cell in the final cell configuration (at time t = 1010 time steps) was born a time t = 6 time steps, and the youngest cells were born at time t = 1009 time steps. The average age of the all the cells in the final cell configuration was 538.5 time steps.
Preliminary work strongly suggests that the control shown is not time optimal. However, in improving the duration we would create a control with more switchings and while the geometric distance would still be 12 the boundary of the final mass would not look as smooth as in this simulation. Since the user set the fractone locations, the control was not guided by any algorithm beyond judging what direction had the most potential for growth without switching fractones at any given time that the model was paused. This control is not unique, either. Different controls could have produced the same final target configuration (although the timing of the growths may have been different). Fractone positions could have been changed at every single time step to produce a final configuration that was exactly the same as the target, however this would be contrary to our hope to find a control that only exerted itself a small amount.
Moving forward with simulations, we have several goals. First is to model actual biological systems using experimental data. Second, we would like to re-introduce the neutral growth, and examine what techniques and strategies are required in order to maintain reasonable accuracy to the target configuration. As mentioned, the use of negative fractones could reduce the impact of the random neutral growth. And third, we will consider optimization of the growth, finding strategies to minimize time or force fractone changes or some other determined cost. 6. Conclusion & future work. We have presented a new control hybrid model of morphogenesis and demonstrated many of its key features through a numerical version. We have framed the model as a controlled hybrid system with stochastic elements. The important mechanics of the model, including the diffusion and capture of growth factor, the deformation of the cell mass caused by meninges, the ability to implement real-world experimental data, and the control asserted by fractone placement, have been shown. We have demonstrated the effect of several parameter value changes in order to develop a better intuition of the behavior of the model as parameters change. This is in preparation for the next phase of the work when we must tune the parameters around actual biological data. We have also examined the question of controllability on the non-stochastic model and offered some insights into the optimization question. While much has been done in preparation, the biological extent of a fractone's ability to control growth will not be known until further laboratory experimentation is completed. Once experimentally-determined three-dimensional maps of developing rat embryos are available, we can begin testing the fractone hypothesis by initializing the model using a map and comparing the results to maps of embryos which are 24 hour older. This process will also allow us to refine the parameter values we have thus far assumed. This includes the growth factor production functions, the timing of the growth factor production, the threshold values for the positive and negative growth factors, the delay between mitosis events, the diffusion constant, the absorption constants, and the initial growth factor distribution. After comparison to biological images, we hope also to compare our results to existing models of cellular growth and development to measure whether the incorporation of fractones generates any more accurate or realistic growth patterns.
The MATLAB code for the model has also not been optimized. We would like to examine larger biological structures and reduce the size of the time step to produce better results, however with the current code this is excessively time-consuming. Better hardware and more efficient coding would help to deliver effective results. The current approach to numerically approximating diffusion shall also be reconsidered to see if a more realistic method can be found.
Finally, determining optimal strategies for growth and finding optimal controls given any pair of initial and final sets of cell bodies is a very interesting problem. While naive methods may lead to efficient executions, they are fairly slow as they require repeated trial and error. Further study of the system may reveal better insights leading to more refined methods and results, as well as algorithms to automate the selection of good control of fractone positions (as compared to an intelligent user sitting down and setting the control as was done in the rabbit example). Understanding of the optimal methods in the model may lead to better understanding of the biological system and grant better insight into how all lifeforms takes on their particular form. Our model focuses on fractone-directed morphogenesis, where fractones are responsible for accelerated and directed growth of cells. A common question posed then is "What directs the fractones?" While this is currently unknown, it is our hope that understanding the questions of controllability and optimization may lead us in the right direction to answer this question.