1. Introduction

Despite recent progress in tissue engineering, it is still difficult to understand the mechanisms underlying tissue growth and maintenance. To progress in regenerative medicine, it is important to understand how cell proliferation, differentiation, migration and apoptosis are regulated and ultimately lead to the formation of a viable tissue.

Due to the implementation of a Virtual Lab ^[1], it is now possible to formulate and test many hypotheses in biology, particularly on tissue growth.

In a mathematical point of view, a tissue is a set of cells. Therefore, the state of the system is a set of set of cells. Tissue growth can be defined as the set of processes that causes a tissue to develop and converge towards its final structuration.

---

2. State of art

Many computational models exist in this field. Different models are designed according to the molecular, cellular or multicellular level of investigation. The models also differ for prokaryotes, animals and plants. For example, the tensegrity models are used to determine how individual cells change their shape. These models are based on biomechanical forces between cells and the extracellular matrix. The tensegrity models show that elongation of cells may result from adherence to the extracellular matrix ^[2]. However, cell differentiation usually occurs prior to the acquisition of a cell shape.

Artificial regulatory networks have also been used to model morphogenesis. Regulatory networks are defined as a series of regulatory genes and structural genes. They follow a network of rules that determine the evolution of the system. The regulatory networks are intended to generate a simple specific pattern ^[3]. However, it is difficult to explain how an organism emerges from an orderly sequence of gene expression ^[4].

In 1952, A.M. Turing proposed a reaction-diffusion model in which a system of chemical substances, called morphogens, react together and diffuse through a tissue ^[5]. In this model, patterns or structures emerge from instability of the homogeneous equilibrium. Random disturbances trigger the instability. This model produces structures such as spots and stripes but it does not provide any information related to the process of morphogenesis. The model based on random disturbances cannot explain in a deterministic way how cell divide, migrate and differentiate from the first cell to the whole organism.

Another model highlights the role of fractones ^[6] as captors or activators of growth factors which allow a much finer control than simple diffusion of growth factors. This could explain the polarisation and the directional divisions of the cells. The fractones are extra-cellular matrix structures that bind, promote and regulate growth factors. It seems that they determine how and when cells divide, migrate and differentiate from the first cell to the whole organism.

The model of fractones is used to describe the growth of shapes ^[7,8] and a hybrid automaton is implemented. Questions remains to understand morphogenesis, why, how and where fractones are placed in the extra cellular matrix, how do they control in space and time cell dynamic to maintain the shape.

Another approach commonly used for modeling multicellular systems is to integrate the dynamics interactions between different spatial and/or temporal scales. This approach introduces some complexities into models that limit their effectiveness. Thus, we must focus on the well, which determines both causalities and downgrades. In other words, as a first step to a better understanding, observations must be restricted to a single spatial and temporal scale of the biological organization.

Thanks to advances in microscopy and imaging, very detailed data on components and structures of living organisms are now available. For example, Melani et al. achieved a tracking of cell nuclei and the identification of cell divisions in live zebra fish embryos using 3D+time images acquired by confocal laser scanning microscopy ^[9].

On the basis of these observations, we formulate a set of principles to propose a model. First, the noted geometrical division allows us to adopt a discrete model in time and space to study the structures appearing in the early stages of morphogenesis. For now, we restrict to a 2D model. We hypothesize that, based on spatial choices of cells, we can define a morphological dynamic of tissues emergence and growth at the early stage of the embryo during which mechanical influences are insignificant. Furthermore, we assume that this dynamic is the fundamental principle of morphogenesis and is therefore able to describe all evolutions of a tissue, both those that modify it, and those that maintain it.

The problem addressed in this article is to understand the cellular mechanisms that allow cells without central control to maintain a shape (tissue, organ, ...) despite renewal of the cells (apoptosis, mitosis) that compose it. This is what we call dynamic equilibrium.

The research work presented in this paper is mainly motivated by the following question : how do geometric cell divisions allow creation and maintainance of a specific shape? We address this question using the mathematical framework of viability theory and mutational analysis ^[10,11,12]. We will use a mathematical model of multicellular developmental design based on morphological analysis ^[13].

Viability is a very important feature of dynamical systems under state constraints whose initial value problem does not ensure uniqueness of solutions. If uniqueness of solutions is not required, then at each time, all solutions have their values in the fixed constrained set (invariant set), or at least, one solution with this property exists (viable set).

Viability theorem is adapted to morphological inclusions ^[14]. Using this mathematical framework, we investigate the set of 2D shapes with 16 cells ^[15] for which we know the states of the cells and the controls they use to make this shape. We compute the viability set at a certain time horizon with an energy constraint.

This article is organized as follows.

First, we present the mathematical morphological dynamic of a multicellular system, which offers an approach especially well suited for tissue evolution. The dynamics model allows us to generate all possible phenotypes of the tissue after 4 division cycles and their underlying genotypes.

Then, because we are focusing on tissue growth, cellular exchange is taken into account through a coupling of an energy diffusion model with the tissue model. Thereby, cells are endowed with the ability to consume, lose, and transfer energy from and to their environment. This model of cellular growth is able to describe tissue evolutions and shape maintenance with resources constraints.

The hybrid automaton model examines the continuous and discrete aspects of the system and how they interface : differential equations manage the energy of the system and morphological equations govern its growth. Coupling these two models implies a coevolution in the system. On one hand, the global energy fed in the system depends on the size of the tissue. On the other hand, the local morphological dynamics of cells that locally transforms the tissue depends on their available energy reserve.

We present here the algorithm of the viability set of shapes with this dynamic under resources constraints at different time horizons and we discuss the results. We also present another extended viability algorithm with all the possible genes completions. These algorithms are implemented in a hybrid automaton using the C++ Boost library for the representation of data and for parallel computation. Finally, we propose an evolutive algorithm passing from one stable shape to another one by changing the control regulation.

---

3. Mathematical morphological dynamics

In many biological articles, we can read that cells create a precise spindle orientation during the first steps of morphogenesis ^[16,17]. Based on this in vivo observation, we propose a mathematical model based on dynamical systems. In a previous work, we introduced the mathematical model of cellular division using mutational and morphological equations ^[11,18,12]. These morphological equations were the first step to comprehend morphogenesis as a morphological dynamic converging to a target (the final shape). This is a cellular growth process where each cell is autonomous and regulates itself with epigenetic feedback from other cells. The shape is the target of a multivalued dynamical system ^[13].

We formalize the dynamics of cell development, i.e. the series of choices or decisions (migration, division, apoptosis) made at each step by cells. We then code this series of choices for all cells to understand how cells with the same code can differentiate and change their behavior depending on the environment. With time this evolution leads to a strong, stable shape.

For each cell, we select only one set of 8 possible actions :

$\{\uparrow, \downarrow, \leftarrow, \rightarrow, \nearrow, \swarrow, 0, \emptyset \}.$

or

$\{1, 2, 3, 4, 5, 6, 7, 8 \}.$

1. $\uparrow$ is the migration moving up from the cell.

2. $\downarrow$ is the migration moving down.

3. $\leftarrow$ is the migration moving to the left of the cell.

4. $\rightarrow$ is the migration moving to the right of the cell.

5. $\nearrow$ is the migration moving to the back of the cell.

6. $\swarrow$ is the migration moving to the front of the cell.

7. $0$ is the quiescence.

8. $\emptyset$ is the apoptosis.

We have to order these actions (128 different ordering sequences in 2D and 3912 in 3D). These permutations of polarity division differentiate the cell controls. In turn, the different cell controls lead to different shapes. At each time the two daughter cells can change the sequence of actions and then differentiate.

This mathematical formalization allows us to reach any shape and provides a general framework to code all possible dynamics. The difficulty is to correlate the code and the shape, ultimately the genotype and phenotype.

---

3.1. Notations

1. $\mathcal{K} \subset \mathcal{P}(X)$^* denotes the morphological environment^† $(X=\mathbb{R}^{3})$, space of "cells" $x \in X \cup \emptyset$ characterized only by their position for the living cell or by death. The shape consists of "tissues" $K \in \mathcal{P}(X)$, $L \in \mathcal{P}(X)$, etc... which are themselves subsets of cells.

^*Supplied with the structure of max-plus algebra for the operation $\cup$ et $+$ (where $K+\emptyset:=\emptyset$).

^†For instance, $\mathcal{K}:= \left\{ K \subset M \right\}$ is the family of subsets contained in a given subset (confinment) $M$.

2. The set of eight (genetic) actions

$\mathcal{A} := \{\uparrow, \downarrow, \leftarrow, \rightarrow, \nearrow, \swarrow, 0, \emptyset \}$

made of the six "geometric directions" plus quiescence and death, are used to describe the transitions between two successive stages.

(a) (geometric) transitions $x \mapsto x+d$, where $d \in \mathcal{A}$ (actual action).

(b) stationarity, $x \mapsto x+0=x$ (no action).

(c) apoptosis; $x \mapsto x+ \emptyset=\emptyset$ (cancellation or suicide).

remark : controls $\emptyset$ or $0$ are always possible. Therefore, the cell are never inactive.

---

3.2. Cellular division

This injunction $(d^{\leftthreetimes}, d^{\rightthreetimes})$ is described by the genetic inclusion

$x \leadsto \left\{ x+ d ^{\leftthreetimes}, x+ d^{\rightthreetimes}\right\}$

(1)

where the mother cell $x$

• first migrates from $x$ to $x+ d ^{\leftthreetimes}$ using the migration action $d ^{\leftthreetimes} \in \mathcal{A}$ at a new position (including $x$ or $\emptyset$);

• and, next, divides giving birth to a daughter cell at position $x+ d^{\rightthreetimes}$ using the $d^{\leftthreetimes} \in \mathcal{A} \setminus \{0\}$.

and produces two daughter cell pair

$\left\{ x+ d ^{\leftthreetimes} , x+ d^{\rightthreetimes}\right\}$

Hence the elementary transitions are described by

1. sterile migration by taking $d^{\leftthreetimes} \in \mathcal{A}$ and $d^{\rightthreetimes}=\emptyset$;

2. stationary division by taking $d^{\leftthreetimes}:=0$ and $d^{\rightthreetimes} \in \mathcal{A}$;

3. migrating division by taking $d^{\leftthreetimes} \in \mathcal{A} \setminus \{0\}$ and $d^{\rightthreetimes} \in \mathcal{A} \setminus \{0\}$.

---

3.3. Gene and status expression

To define the genetic process, we will introduce a set of permutations $\sigma$ of the set $\mathcal{A}$.

A genetic process is an ordered sequence of actions

$d_{\sigma}:= \{d_{\sigma (1)}, \ldots, d_{\sigma (8)}\} \in \mathcal{A}$

We define by $\mathcal{G}$ the subset of genetic processes of actions in $\mathcal{A}$ which is itself the subset of eight permutations.

The order of actions will highly impact the shape at the early stage of growth. This will ultimately determine the final shape of the tissue.

Operating a genetic process under a given criterion, either for migration or for division, means that the process scans successively $x+d_{\sigma (1)}, \ldots, x+d_{\sigma (8)}$ up to a possible transition.

All genetic process will lead to the empty set $\emptyset$ (apoptosis) or $0$ (quiescence). Therefore, the scanning operation will stops before the eight genetic actions. Consequently, there are only 128 different ordering sequences in 2D and 3912 in 3D.

---

3.4. The Double Time Scale

In order to define the morphological evolutionary mechanism $K_{n} \mapsto K_{n+1} := \Phi (n, K_{n})$, we will distinguish :

• the calendar (or algorithmic) time $n \geq 0$,

• the process time $j=1, \ldots, j_{m}$, at each time $n$, for intermediate operations mapping $K_{n}$ to $K_{n+1}$. At each process time, the daughter cells cannot occupy a location already taken by another cell.

A genetic regulon is defined as a map^‡ that associates each triple $(n, L, x)$ to the pair $(G^{\leftthreetimes}(n, L, x), G^{\rightthreetimes}(n, L, x)) \in \mathcal{G} \times \mathcal{G}$. The genetic processes satisfies the following non-overlapping properties:

$\left\{ \begin{array}{ll} \forall \; x \in L, \\ i \;\mbox{ is the first stage in the}\; genetic\; process\; \mbox{such that}\; x +G^{\leftthreetimes} (n, L, x)(i) \; \in \; \{x\} \cup \complement L\\ j_{i} \;\mbox{is the first stage in the} \; genetic\; process \; \mbox{such that}\; x +G^{\rightthreetimes} (n, L, x) (j_{i}) \; \in \;\{x\}\setminus \{x+G^{\leftthreetimes} (n, L, x)(i)\}\cup \complement L \end{array} \right.$

(2)

^‡A single-valued map for the time, since no other (biological) parameters are involved at the stage of this study.

The first property determines the migration of the first cell such that $x +G^{\leftthreetimes} (n, L, x)(i) \in \complement L$.

The second property results into the division and the birth of the two daughter cells. The first cell is at the position of its mother and the second cell will occupy a free space $G^{\leftthreetimes}(n, L, x)(j_{i})\ne G^{\leftthreetimes} (n, L, x)(i)$).

1. Apoptosis occurs whenever $G^{\leftthreetimes} (n, L, x)(i)=\emptyset$. In this case, cells will not divide;

2. Cell migration will not occur whenever $G^{\leftthreetimes} (n, L, x)(i)=0$. In this case, $x +G^{\rightthreetimes} (n, L, x) (j_{i}) \in \complement L$;

3. Cell division will not occur whenever $G^{\rightthreetimes} (n, L, x) (i)=\emptyset$.

---

3.5. Local Morphological Dynamics

The genetic regulons $(G^{\leftthreetimes}, G^{\rightthreetimes})$ are either located in the first cell or constructed to obtain a viable organism. The genetic regulons will determine the following controlled morphological inclusion :

• First, the map $H$ defined by

$H(G^{\leftthreetimes}, G^{\rightthreetimes})(n, L, x) \; := \; ( x +G^{\leftthreetimes} (n, L, x)(i), x +G^{\rightthreetimes} (n, L, x) (j_{i})).$

(3)

It associates with any pair $(G^{\leftthreetimes}, G^{\rightthreetimes})$ of genetic processes to the daughter cells.

• Second, the transition transforms the subset $L$ at time $n$ for the cell $x \in L$ to the subset

$\varphi (n, L, x; G^{\leftthreetimes}, G^{\rightthreetimes}) \; := \; L \bigcup\limits_{}^{} \{G^{\leftthreetimes}(n, L, x), G^{\rightthreetimes}(n, L, x)\}$

---

3.6. Global Morphological Dynamics

We assume that a subset $K_{n} \in \mathcal{P}(x)$ is coded by an ordered list $(x_{1}, \ldots, x_{p_{K_{n}}})$.

Hence we construct $\Phi$ in the following way.

Given $K_{n}$, we define the sequence :

1. $K_{n}(x_{1}) \; := \; \varphi(n, K_{n}, x_{1}; G^{\leftthreetimes}, G^{\rightthreetimes})$

2. $\forall \; p = 2, \ldots, p_{K_{n}}$,

$K_{n}(x_{1}, \ldots, x_{p}) \; := \; \varphi(n, K_{n}(x_{1}, \ldots, x_{p-1}), x_{p}; G^{\leftthreetimes}, G^{\rightthreetimes});$

(4)

3. We set $\displaystyle{\Phi (n, K_{n}) \; := \; \varphi (K_{n}(x_{1}, \ldots, x_{p_{K_{}n}}), G^{\leftthreetimes}, G^{\rightthreetimes})}$.

This mathematical formalization is used in a virtual reality simulator. The cells have the same virtual genetic material which is a list of $(G^{\leftthreetimes}, G^{\rightthreetimes})$.

---

3.7. Reachable sets of phenotypes

This formalization of cell growth is based on autonomous cells that perceive the microenvironment. The cells use the control part of their genetic heritage to maintain their viability or to continue their growth after a punctuated equilibrium. This model was developed in 3D using a multi-agent based model ^[1]. The instructions are contained in a genetic program in the first cell. This is the only instruction given to the cells.

A big challenge is to determine which shapes are biologically viable.

To address this question, a 2D-cellular automaton is implemented ^[15]. A systematic search of reachable shapes is done for all possible dynamics. We have captured all the reachable shapes from a single cell to 16 cells (after 4 division cycles). Once a new shape is created, we identify its creation process, the phenotype, and the genotype.

Starting from this set of shapes with known genotype, we will look at the shapes that remain viable at a given time horizon. This defines the viability kernel. We will study the dynamic of cells, the available resources and the stability of cells that construct and maintain the viability of the biological structure.

---

4. Prediction tool based on the links between cell fate and their access to energy

To model growing tissues, we introduce energy management and exchange between cells and their environment. This allows us to define conditions under which cellular mechanisms such as quiescence, mitosis, and elective cell death are held (see Fig. 1).

To define $Init$, the initial state of the tissue, the same amount of energy is given to each cell. Each cell has access to a quantity equivalent to $E_ {ma}$. In the absence of interactions, the concentration of energy accessible is the same for all cells.

These conditions refer to three available energy thresholds that have to be initialized: minimum $E_{mi}$, medium $E_{me}$, and maximum $E_{ma}$.

We define $E=n(x)$ as the available energy of cell x at the beginning of a given cell division $n$. Prior to reading its genetic process $G_x$, $x$ compares its current available energy $E(x)$ to these thresholds, then the cell behaves appropriately by choosing the corresponding mechanism.

The correspondences between threshold and mechanism are fully described in ^[15].

The dynamics of the hybrid system $F$ is given by the variation of $E_n (x)$ and by the morphological dynamics of $K_n$. This model is implemented in the hybrid automaton.

---

5. Viability kernel algorithm

Viability theory deals with evolution equations and proposes concepts and methods to control a dynamical system in order to maintain it in a set of viability constraints. There are many applications in the field of ecology, economy or robotics. For example, when a system leaves an area of its state space and then dies or deteriorates.

The growth of a tissue under certain conditions such as energy constraint that the cells receive can be viewed as a viability problem with multivaluate equations.

Viability kernel algorithms have been developed to approach the viability kernel. For example, in ^[19], the authors use an approximation of the viability kernel by a discrete sequence of viability kernel. Simirarly, other authors have developed an algorithm that computes a progressive approximation of the viability kernel using a classification method based on support vector machines ^[20]. In another study, a numerical approximation of the viability kernel was computed using the tools of dynamic programming ^[21]. A viability kernel algorithm is used to find the viable controls at each step to avoid an exponential increase in computation time as a function of the size of the control space. Even for a continuous system, all these algorithms for the approximation of the viability kernel use a grid to define a discrete system and calculate the exact kernel of the system.

To our knowledge, there is no algorithm to calculate the viability kernel applied to the biological shapes.

Since a biological shape is viable for a given time, we will study the viability kernel up to a time horizon. The viability kernel is the set of all states of a system from which at least one control function ensures compliance with constraints up to this time horizon.

We will therefore calculate the viability kernel of all tissues reached after 4 cycles of division. In other words, we determine the subset of viable tissues that maintain their shape and respect ressource constraints up to a certain number of cycles without acquisition of new gene.

---

5.1. Variables and Parameters

To calculate the viability kernel for this problem, we define :

• The state of the system, describing the phenotype.

    – The set of shapes after 4 divisions cycles (1029 tissues with 16 cells) ^[22].

• The time horizon: this is the number of cell cycles for which the tissue remains viable.

• The controls of the dynamics of the shapes, i.e. the genotypes.

    – The dynamics of tissue growth is ensured by cell's genome.

    – This genome must ensure a viable development up to the time horizon without acquisition of new genes.

• Constraints express the conditions under which the state of a system is considered viable.

    – The constraints on tissues refers to the level of energy above a threshold for its maintenance.

    – Example of constraint: the level of energy that maintain 50% of the cells.

---

5.2. Dynamics

Developed software using classical viability kernel algorithms exist for controlled differential equations. For morphogenesis, we implement a new viability kernel algorithms for controlled mutational equations.

For the computation of the viability kernel of a tissue, a state is a set of points representing the cells. The evolution of the tissue depends on cell migration, proliferation or apoptosis, and this for each cell of the tissue, at all stages of the tissue formation up to a certain time horizon. At each cycle $t$, we start by the set of tissues in the current viability kernel ($Viab_t$) and we compute the $t+1$ kernel ($Viab_ {t +1}$).

There are two problems of complexity in this computation :

• the algorithmic complexity.

• the visualization of the tissues.

We choose a persistent storage of the states, this allow a calculation from a previous kernel. Therefore, informations associated with the tissues are read and modified directly from files.

The time step in the calculation of the viability kernel corresponds to the time required to scan throughout the tissue space. We scan all the possible behavior for each cell of each tissue of the set.

Viability algorithm is less time-and memory-consuming than the calculation of reachable sets. The latter apply for all cells of all tissues of the set, all possible actions. For the reachable set, we have an exponential growth of the number of tissues as a function of the number of cells and for the viability algorithm, it is a linear growth as a function of the number of cycles.

---

5.3. Algorithm

We describe the steps of the algorithm in figure 2. A vector identifies the names of the files that define the tissues of the state space. Thus, as long as the number of cycles is not reached, we will :

• Scan the vector of tissues in order to load the file.

• Simulate the tissue and if its energy level respects the constraints, add the tissue to the vector containing viability kernel. If the constraints are not respected, the tissue is removed from the viability kernel and another tissue is chosen.

• If we have finished the scanning of the state space, it means that all the tissues of a cycle of division have evolved. We transfer the names of the viable tissues for this cycle in the vector and we increment the number of cycles.

• We stop the calculation of the viability kernel in the 3 following cases :

    1. We have reached the number of cycles defining the time horizon. The viability kernel is the kernel computed at the current division cycle.

    2. The viability kernel is invariant for numerous cell cycles.

    3. The computed viability kernel is empty. In this case, the viability kernel is the last set computed at the previous cell cycle.

1. Initialization     $Viab_0 = \{L_0\}, t=0$ 2. Loops     $Viab_{t+1} = \{L_t \in Viab_t~ | ~E(L_{t+1}) > E^* \text{ Maintenance energy}\}$,     $t=t+1$ 3. Stop condition       • $t = t^*$ (The time horizon) or       • $Viab_t = \emptyset$

---

6. Results of viability kernel algorithm

We have computed the viability kernel at different time horizons :

• At 10 cycles, the volume of the viability kernel is 1013 (see figure 3),

• At 100 cycles, the volume of the viability kernel is 225 (see figure 4),

• At 200 cycles, the volume of the viability kernel is 220 (see figure 5),

• At 1.000 cycles, the volume of the viability kernel is 70 (see figure 6),

• At 10.000 cycles, the volume of the viability kernel is 5 (see figure 7).

The 5 tissues of $Viab_{10.000}$ (see figure 8) are those that can maintain a level of energy for $10.000$ cell cycles.

There is a value of maximum energy parameter for each of the 5 tissues of $Viab_{10.000}$ (see table 1). Beyond this value, the tissue converges to a stable form (see figure 9) and indefinitely remains in the viability kernel throughout all time horizon.

In the space of shapes, the stable shapes are limit cycles of this dynamical system. In our case, two different shapes $(427,445)$ with 16 cells converge to the same periodic limit cycle. Furthermore, there are some values of $E_{ma}$ below which the tissue will not maintain its shape. These values change the behavior of the dynamical system like bifurcation points.

In a previous article, we performed preliminary simulations ^[22] to understand the effect of radiotherapy on tumor. The efficiency of the therapy depends on the tumor shape. Irradiation applied to stable shapes might be detrimental. This new approach compelled us to classify stable shapes at a long time horizon and to study the link between the level of energy and the number of cycles.

When $E_ {ma} < 0.0502$, $Viab_ {33} = \emptyset$. Similarly, a value of $E_ {mi}$ warrants a convergent evolution that indefinitely remains in the viability kernel. Indeed, if $E_ {mi}> 0.020885$, $Viab_5 = \emptyset$ (see table 2).

Beyond a level of energy for each tissue, the evolution leaves the viability kernel. On the other hand, the tissue converges and undergoes a perpetual evolution when this level of constraint is reached (see table 3).

To study the convergent and perpetual evolution of these tissues, we have traced the variations of their shape through 10.000 cycles. When these tissues have reached their convergence cycle, they evolve through specific forms with a periodic evolution.

Some animals use this periodic mechanism. For example, turritopsis dohrnii repetitively turns into cycle of adult and polyp throughout his live and never die.

The periodic changes in the tissues are shown in figures 10 to 14.

The periodic evolution is a particularity of the 5 tissues revealed by $Viab_{10.000}$. Any other tissue in the set does not display this characteristic.

The figure 15 shows the variation of the of the tissue shape 249 through the cycles without periodicity.

---

7. Extended viability algorithm

In the previous computation of the viability kernel, the evolution of the states is associated with the controls used to reach the shape. If the resulting dynamics do not respect the constraints, then the state leaves the viability kernel. However, the most general case of the computation of the viability kernel is to consider the set of all possible evolutions. This is equivalent to apply the set of possible controls to the state.

When the genome acquired in the first division stage (first 4 division cycles) is incomplete and the viability constraints are no longer respected (crisis of viability), the genome of the cells can be enriched without the acquisition of new genes. The last gene is completed with all the possible combinations of the missing genetic actions. When the tissue exit its viability constraints, the last gene is completed with all the possible combinations of the missing genetic actions.The objective keep the tissue in the viability kernel.

The generalization of the viability kernel computation consists in testing all the possible controls in order to find a viable dynamics for the tissue. If one of the control is viable, the tissue with this completed genome is in the viability kernel. On the other hand, if there is no genetic complement that allows the tissue to stay in the viability kernel, the tissue is removed from the set of viable tissues.

The algorithm Fig. 16 is as follows:

• We will retrieve the last genome of the gene

• Then, we identify all the missing genetic actions

• In order to have exhaustively all the possible genetic completions, we make prints with order and without repetition in the list of missing actions. The number of draws is determined by the number of missing actions. For example, if the complement contains only the genetic action 4. Then the missing actions are 1, 2 and 3.

• So, we will proceed to three draws with order and without repetition, which gives 15 genetic complements (see figure 17).

• Thus, by making a tissue evolve, if the dynamics is not viable, we take up the dynamics by modifying the genome with the addition of a genetic complement. As long as it is not viable, we continue to test other genetic complements, until the dynamics are viable or we have no more complements to be tested.

• As soon as there is a viable dynamic, we consider that the state belongs to the viability kernel. On the other hand, if all the dynamics leave the constraints space then the state does not belong to the viability kernel.

This generalization naturally increase the volume of the viability kernel. In addition, we have for this simulation an extra cost generated by the search for a viable evolution by testing all the genetic complements.

The new viability kernel is calculated at time horizons 10,100,1000 and 10.000 cycles. The results showed an enrichment of the viability kernel with a larger volume because at each time horizon, for each state, all the possible dynamics are attempted before excluding the state of the viability kernel. For some states in some time horizons, it has resulted a viable dynamic.

• At 10 cycles, the volume of the extended viability kernel is 1025 (see figure 18),

• At 100 cycles, the volume of the extended viability kernel is 402 (see figure 19),

• At 1.000 cycles, the volume of the extended viability kernel is 120 (see figure 20).

• At 10.000 cycles, the viability kernel does not change with respect to the previous version of the algorithm.

The viability kernel of stable shapes did not grow. The number of stable shapes will grow only when a new control is introduced.

We have compared the volumes of the viability kernel between a calculation with a single possible dynamic for each state and a calculation with several dynamics (see figure 21).

---

8. Discussion

The volume of the viability kernel is influenced by the control space and other factors. Raising the constraint requirements to 90 % caused the volume of $Viab_ {10}$ to decrease to 1010. 15 tissues no longer satisfying the new constraints were excluded from the kernel. When the constraints were reduced to 10, the volume of $Viab_ {100}$ increased to 563. Then, 161 tissues satisfied the new constraints. Moreover, by varying the thresholds of energy, we can also observe differences in the volumes of viability kernels.

We can already observe the influence of the variation of the thresholds of energy on the morphologies of the tissues. In figure 22, we show the tissue shapes 69 with 10 division cycles for different values of $E_ {ma}$.

For the calculation of reachable sets, we encountered the difficulty associated with the representation of the reachable shapes and the tree walk. The constraints strongly decreases the number of shapes and open new perspectives for the application of this theory to biological shapes evolution.

The viable set revealed that there are tissues that can remain viable ad-vitam aeternam regardless of the time horizon. This concept of dynamic equilibrium is crucial to understand the evolution of shapes. How an organism evolve if it hasn't reached equilibrium? Starting from these shapes, multicellular organisms react to external stress and change their regulations in order to reach another viable equilibrium. This viable mechanism evolves shapes from simple viable structures into more complicated viable structures.

The stable shapes are viable shapes for a given biological constraint. Depending on the constraints, viables shapes can change from one "stable" shape to another. In the space of shapes, this equilibrium of shapes are like periodic limit cycles for a given dynamic. With the adequate biological constraints we could restrict our set of evolving shapes from these stable shapes and then at each step compute the new viable set of shapes with new set of controls.

---

Conflict of Interest

All authors declare no conflicts of interest in this paper.

---