Forest structure in epigenetic landscapes
Abstract
Morphogenesis is the biological process that causes the emergence and changes of patterns (tissues and organs) in living organisms. It is a robust, selforganising mechanism, governed by Genetic Regulatory Networks (GRN), that hasn’t been thoroughly understood. In this work we propose Epigenetic Forests as a tool to study morphogenesis and to extract valuable information from GRN. Our method unfolds the richness and structure within the GRN.
As a case study, we analyze the GRN during cell fate determination during the early stages of development of the flower Arabidopsis thaliana and its spatial dynamics. By using a genetic algorithm we optimize cell differentiation in our model and correctly recover the architecture of the flower.
keywords:
Biomathematics, Gene regulatory networks, Epigenetic landscapes, Morphogenesis model1 Introduction
Genetic Regulatory Networks (GRN) govern morphogenesis, the process responsible for producing the complex shape of full grown organisms from a few cells (KellerMorph (), EttensohnMorph ()) and every process of life including metabolism, the cell cycle and signal transduction. Genetic Regulatory Networks (GRN) are a composed of a collection of genes that interact between each other and with external factors, forming a complex network (DavidsonGRN ()).
In recent years GRN have gained a lot of attention and great advances have been made (KarlebachGRN (), EmmertGRN ()). From cancer identification (IyerCancer ()),(ChangCancer ()) and diabetes (JesminDiabetes ()) to the formation of animals’ bodies (DavidsonAnimal ()) GRN are being used to model and understand diverse biological processes. However these complex networks haven’t been thoroughly understood and there is a lot of experimental data that needs to be analyzed using computational and mathematical techniques.
In this work, we propose Genetic Regulatory Trees, as a way to study Genetic Regulatory Networks. By modeling GRN as a discrete dynamical system (boolean automata) and analyzing all its possible states, we find that they have a forest structure, where each tree corresponds to a fixed state of the automata. This forest models Waddington’s epigenetic landscape of the network (a developmental model that illustrates the mechanics of cell fate differentiation see section (2.2)). The use of this technique to model GRN and its epigenetic landscape is a novel approach.
As a case study we analyze the GRN of the flower Arabidopsis thaliana during cell fate determination. In (AlvarezFM ()), using experimental data, we obtained the gene regulatory network (GRN) of the flower Arabidopsis thaliana that determines the fate of floral organ cells. The network models the interactions between genes responsible for cell differentiation: each node in the automata represents a certain gene, the edges are the interactions between them and the logical rules that govern the automata are based on detailed experimental data.
Genetic Regulatory Trees are defined and built as a directed graphs that turn out to have a directed rooted tree structure. Each node in a tree represents a given state of the Boolean automata (GRN) and every possible state is considered. The edges link states that follow each other in the dynamical system and the root of each tree represents a steady state of the automata, so that there will be as many trees as steady states. Once the trees are built, we define chains of cell types that will traverse the trees in a certain order. Each link in the chain will represent a group of undifferentiated cells (of the meristem) with the same genetic configuration, and neighbouring links will be genetically similar to each other. Biologically we can think of it as a radial longitudinal section of the meristem. We define an energy measurement of a chain (the energy it requires for all of its links to differentiate). Considering nature is efficient, we use Genetic Algorithms to minimize this energy and expect to recover accurately the spatial configuration of the flower during cell fate determination. The set of trees represent the epigenetic landscape that models how the different environmental and genetic forces affect cellular differentiation.
This method is general enough to be used in any GRN as a tool to retrieve very interesting information about them, like the robustness of cell differentiation into each flower organ (trees), the importance of a specific state (a node in the tree), and the underlying dynamics of certain processes, like cellular differentiation, as studied in this work.
This paper is organized in the following manner: The first section is this introduction; on section two we present biological background information and we detail the discrete dynamical system (the Boolean automata) of the GRN we will analyze. On the third section we present the details of our model and we construct our Genetic Regulatory Trees taking into account experimental data (from the discrete dynamical system). On section 4 we present the details of the genetic algorithm, used to minimize the energy of the paths, and the obtained results. Finally on the last section we discuss our concluding remarks.
2 Background information
It this section we provide the biological information needed to understand our work and describe the Genetic Regulatory Network model we are using.
2.1 Biological background
The flower organs of all species of Angiosperms (approx. 250,000) are organized in four concentric rings (whorls), which are, from the outer rim to the center: sepals, petals, stamens and carpels (fig (1)). The only known exception to this configuration is the one observed in the flower Lacandonia schismatica where the position of its stamens and carpels is inverted.
We work with the flower of the plant Arabidopsis thaliana. This plant was the first one whose complete genome was sequenced and has been extensively studied (ArabidopsisGen ()), (AlvarezArabidopsis ()).
In (AlvarezFM ()), the gene regulatory network (GRN) of the flower Arabidopsis thaliana was obtained. It is based on experimental data and it’s defined as a Boolean automata which determines the fate of floral organ cells.
Using this model, we define and construct our Genetic Regulatory Trees as we will see in the following section.
2.2 Epigenetic landscape
Epigenetic landscapes, originally proposed by Waddington in 1975 (WaddingtonEL ()), are developmental models that illustrate the mechanics of cell fate differentiation. These models use as a metaphor a mass in a potential field with a certain number of basins of attraction and paths that lead to each one of them. We propose a way to model the epigenetic landscape associated with a certain GRN using Genetic Regulatory Trees. As a case study we use the GRN of the flower Arabidopsis thaliana during cell fate determination. Each tree will correspond to a different flower organ (sepals, petals, stamens, carpels and the inflorescences). We use the information obtained from the Boolean network, a discrete dynamical system to build them. In the following section we present the details of this network.
2.3 Boolean network
In (AlvarezGRN ()) a discrete dynamical system was used to explore the dynamics of cell fate determination during the early stages of flower development. The system is a Boolean automata, consisting of nodes, each one corresponding to a specific gene ^{1}^{1}1The genes considered are: FUL (Fruitfull), FT (Terminal flower), AP1 (Apetala 1), EMF1 (Embrionic flower 1), LFY (Leafy), AP2 (Apetala 2), WUS (Wuschel), AG (Agamous), TFL1 (Terminal flower 1), P1 (Piscilata 1), SEP (Sepallata), AP3 (Apetala 3) and UFO (Unusual flower organ) AlvarezArabidopsis ()..
Each node has two possible states ( or ) that are updated according to experimentally obtained rules, that correspond to the interaction between genes. These interactions can be expressed as a directed graph as shown in figure (2), where we can see a green edge from node to node if node ’s state affects and a purple edge between node and if both nodes affect each other. The rules that the automata follows (obtained experimentally) can be found in (AlvarezGRN ()).
There are possible initial conditions (binary numbers form to ). The system is iterated, starting from each initial condition. It converges to eight different attractors ^{2}^{2}2For simplicity we don’t take into account two additional attractors that are irrelevant in this work because they correspond to flower organs that are already considered (petals and stamens) and the number of initial conditions that reach them is very small.. Every single initial condition lands in one of these ten states. Each atractor represents one of the main cell types observed during the early stages of flower development (the meristematic cells of the inflorescence and the primordial cells of the flower meristems of sepals, petals, stamens and carpels) EspinosaGRN (). Each equilibrium point has thirteen components (see table(1)).
Floral organ Atractor  

Inflorescence 1  
Inflorescence 2  
Inflorescence 3  
Inflorescence 4  
Sepals  
Petals  
Stamens  
Carpels 
The number of nodes (initial conditions land in each fixed point) in each tree will be the following:
(1)  
(2) 
where the subscripts S, P, E and C correspond to sepals, petals, stamens and carpels respectively and to the inflorescences. (EspinosaGRN (),AlvarezFM ()).
Based on this discrete model, we construct the Genetic Regulatory Trees.
3 Our method
3.1 Genetic Regulatory Trees
Genetic Regulatory Trees model the epigenetic landscape of Genetic Regulatory Networks and are constructed by analyzing the states of the discrete model detailed in the previous section.
The Genetic Regulatory Trees are directed graphs that turn out to have a directed rooted tree structure. To build them we proceed as follows: Each node in the tree will correspond to a specific state of the discrete dynamical system detailed above. Every possible state will be considered (in our case there will be a total of nodes in the trees). There will be a directed edge from node to of the tree if (and only if), following the rules of the automata, the state of node leads to the one corresponding to node (see fig. (3)).
Proceeding this way until every state is analyzed, we obtain an intree (orientated towards the root) for every fixed state in the discrete dynamical system. The root of each tree will be precisely a fix state, that in our case represents a flower organ. The initial points on the dynamical system, that is, the points which are not the iterate of any other point will be the leafs of the trees.
Since our dynamical system has fixed states, we will end up with different trees, which can be seen in the following graphs:
These trees represent the epigenetic landscape of the Genetic Regulatory Network. As we can see, their structure is very interesting.
3.2 Chains of cell types
We are now ready to define chains of cell types. A chain will be formed by links, each of which will represent a group of undifferentiated cells with the same genetic configuration (a possible state of the automata), that is, an array with entries (one for each gene) that can be in one of the two possible states we are considering ( or ). We’ll construct our chains so that the first link will be a node on the tree corresponding to sepals and neighbouring links will be genetically similar to each other.
Biologically, we can think of a chain as a radial longitudinal section of the meristem (recall that flowers are radially symmetrical) (see figure(12)).
The formal definition is the following:
Definition 3.1
(Undifferentiated cells chain) A chain of undifferentiated cell types will consist of an array of links,
where each link is a possible state of the automata, that is, and . The initial one, will be a node in the starting tree (sepals) and each consecutive link will be formed by modifying randomly one (and only one) of the entries of the previous one.
This construction guarantees that neighboring cell types in the chain will be similar to each other. To measure how similar two cells are, we will use Hamming distance:
Definition 3.2 (Hammingdistance)
Given two strings of the same length , the Hamming distance between them is the number of positions at which the corresponding symbols differ (it measures the minimum number of substitutions needed to change one string into the other).
As we can see, Hamming distance between any two neighboring cells in a chain will be one. Note that a cell type (link) in a chain belongs to a given tree (whose root is a steady state of the automata) if it is an initial condition of its steady state. Since in our case, there are nodes in different trees whose distance between them is one, a chain can have links in different trees. We say that a chain traverses a certain number of trees.
Also note that trees are acyclic graphs, so that there will be at most one directed path from any one node to another one, in particular there will (always) be a unique path between from node to the root of the tree.
Definition 3.3 (Specialization path)
The path from a given node in a tree to its root is called the specialization path of and we denote it as .
A chain specializes (differentiates) when each of its links reaches a steady state (following the rules of the automata), that is if belongs to tree , it must traverse its specialization path reaching this way, the root of .
A chain that differentiates correctly (starting with a link at sepals) is the one whose first links specialize to sepals, next links specialize to petals, then stamens and carpels at the end (i.e. it recovers accurately the spatial configuration of the flower during cell fate determination). Recall that an undifferentiated chain is a radial longitudinal section of the meristem, so that the differentiated (or specialized) chain will be a radial longitudinal section of the flower (see figure (13)).
3.3 Energy function
We now define the energy of a specialization path and a function of energy of a chain that measures the energy it needs to differentiate.
Definition 3.4 (Energy of a specialization path)
The energy of the specialization path from , (where is the root of the tree we are on) is given by
Definition 3.5 (Energy of a chain)
The energy of a chain is given by
Considering nature is efficient, we expect the chains that require the least amount of energy to specialize, to be the ones that differentiate correctly. Hence, we will look for the chain of minimal energy. This is a non standard optimization problem that we will solve using a Genetic Algorithm.
4 Genetic Algorithm
We wish to minimize the energy of the chains, that is, we wish to choose the chain that requires the least amount of energy to specialize. For this we design and implement a Genetic Algorithm. Note that a chain with more than links ( being the number of components in each link) is redundant, so it will not be optimal (one with more than links can always be reduced to one with ). We can thus, limit our search to chains with .
4.1 Individuals
An individual of the genetic algorithm together with an initial condition (a starting link) will determine a chain. The starting link of each individual will be a leaf on sepals tree (same value for every individual). From there we will generate the rest of the links by copying the previous one and modifying a randomly chosen component (if it was a it will become a and viceversa). We can represent each individual as
where represents the component that will be modified in the following link, so that the chain generated by will be
where , and in general and (every component but one will be the same between neighbouring links, so that .
We start by randomly generating a population of individuals. will be a fixed number throughout the execution of the algorithm.
4.2 Evaluation and stopping criteria
Having generated a whole population we will evaluate each chain , that is, we will compute its energy as defined in the previous section. Its fitness will be a normalized value, inversely proportional to its energy value in such a way that the lower their energy is, the closer their fitness will be to :
where is a normalization function.
Once the whole population is evaluated we can proceed to apply the genetic algorithm operators (selection, crossover and mutation). These will generate a new population that will again be evaluated. The whole process will stop when we reach a minimum. Since we don’t know the exact value of the minimum, the iterations will stop once the difference between the fitness values of one generation to the next is smaller that a given value .
4.2.1 Selection
Once we have a set of individuals (a population) that have been assigned a fitness value, we will select the individuals that will survive on to the next generation. We use the roulette wheel selection method, a fitness proportionate selection, where each individual has a probability proportional to its fitness value to be chosen. Note that, for the next generation, each individual might be selected more than once, and there might be individuals that are not selected at all. To do this, we start by sorting the individuals by their descending fitness values. The accumulated normalized fitness values are computed (the accumulated fitness value of an individual is the sum of its own fitness value plus the fitness values of all the previous individuals). The accumulated fitness of the last individual should be . A random number is chosen. The selected individual is the first one whose accumulated normalized value is greater than . We repeat this process until individuals are chosen, so that the next generation is complete (the number of individuals in each generation will remain fixed through out all the iterations).
4.2.2 Crossover
We use a two point crossover. For this, we select two individuals and using the selection operator. We will randomly choose two different crossover points . The two individuals will be combined forming two new ones: and . The first part of , up to the crossover point and from to the end, will be copied from , the second part (from until ), will be copied from . Conversely, the first part of (up to ) and from to the end will be copied form and the second part from (see figure (14)). The two new individuals and each define new chains and that will substitute the old ones. We will have a fixed percentage of individuals that will be combined, in our work, we found experimentally that the best value was .
4.2.3 Mutation
The mutation operator is used to maintain genetic diversity from one generation to the next. It allows us to avoid local minimums when working with optimization problems. We will select one individual using the selection operator. We will randomly choose two mutation points . We will replace the value of the states in by a randomly chosen one between and (see figure (15)). This will produce a new individual (with two mutated components):
that defines a new chain that replaces the original one. As with crossover, we select a percentage of individuals that are mutated. We found experimentally to be the best choice.
4.2.4 Elitism
To make sure the we don’t loose the best individual in each generation we apply elitism. We will find the individual with the highest fitness value, and copy it to the next generation with out modifying it. This will guarantee that the best individual in each generation is at least as good as the one in the previous one.
4.2.5 Results
We use the following parameters in our algorithm:

Length of the chain (number of undifferentiated cell types in a chain is

Maximum number of possible iterations allowed:

Crossover probability:

Mutation probability:

Number of individuals per generation:

Elitism:
After running the genetic algorithm an optimum solution is always found in less than iterations (generations). We show an execution where the minimum was found in iterations and the starting point is (a leaf in the sepals tree).
The chain with the lowest energy and the tree that each state belongs to, is shown in table (2).
State  Tree  Specialization Energy 
0101100000000  S  12 
0110110001100  S  1 
0110110001110  P  0 
1110110001110  P  1 
0110110001110  P  0 
0010110001110  P  1 
0010110001111  P  17 
0110110001111  P  0 
0110110101111  T  2 
0100110101111  T  1 
1100110101111  T  0 
1100110101101  T  1 
1100110101100  C  0 
Total energy:  E=36 
where S, P, T, C correspond to Sepals, Petals, Stamens and Carpels respectively.
As we can see, this minimum energy chain traverses the trees in the correct orden, that is, following the flower’s architecture which shows that we are correctly modelling it.
5 Concluding remarks
We are presenting a novel approach to analyze Genetic Regulatory Networks that allows us to unfold the richness and structure within them. In this case we used them to study the Genetic Regulatory Network of Arabidopsis thaliana during cell fate determination and correctly recovered its architecture.
Genetic Regulatory Trees allow us to use tools from graph theory, finite field dynamical systems and finite field arithmetic to analyze the Genetic Regulatory Network behind them and to understand the principles that govern these networks better, which is our next goal. This method is general enough to be used in any GRN.
References
 (1) E. AlvarezBuylla, M. Benítez, A. CorveraPoiré, A. Chaos, S. Folter, A. Gamboa de Buen, A. GarayArroyo, B. GarcíaPonce, F. JaimesMiranda, R. PérezRuiz, A. PiñeyroNelson, and Y. SánchezCorrales, Flower development, The Arabidopsis Book 8 (2010).
 (2) E. AlvarezBuylla, A. Chaos, M. Aldana, M. Benítez, Y. CortesPoza, C. EspinosaSoto, D. Hartasánchez, B. Lotto, D. Malkin, G. EscaleraSantos, and P. PadillaLongoria, Floral morphogenesis: Stochastic explorations of a gene network epigenetic landscape, PLoS ONE 3 (2008), no. 11, e3626.
 (3) Jung H.C. Kim S.H. Lee J.H. Kim J.H. Han S.W. Chang, H.J., Gene regulatory network analysis for triplenegative breast neoplasms by using gene expression data, Journal of Breast Cancer 3 (2017), 240–245.
 (4) Erwin D.H. Davidson, E.H., Gene regulatory networks and the evolution of animal body plans, Science 311 (2006), 796–800.
 (5) Peter I. S. Davidson, E. H., Genomic control process: Development and evolution, Academic Press, 2015.
 (6) M. Benjamin Haibe K.H. Emmert S. F., Matthias, Gene regulatory networks and their applications: understanding biological and medical problems in terms of networks, Frontiers in Cell and Developmental Biology (2014), 2–38.
 (7) C. EspinosaSoto, P. PadillaLongoria, and E. AlvarezBuylla, A gene regulatory network model for cellfate differentiation during arabidopsis thaliana flower development that is robust and recovers experimental gene expression profiles, Plant Cell 16 (2004), 2923 to 2939.
 (8) C. A. Ettensohn, Encoding anatomy: developmental gene regulatory networks and morphogenesis, Genesis (2013), no. 51, 383–409.
 (9) The Arabidopsis Genome Initiative, Analysis of the genome sequence of the flowering plant arabidopsis thaliana, Nature 408 (2000), 796–815.
 (10) Osmanbeyoglu H.U. Leslie C.S. Iyer, A.S., Computational methods to dissect gene regulatory networks in cancer, Current Opinion in Systems Biology 2 (2017), 115–122.
 (11) Jamil H. Hontecillas R. BassaganyaRiera . 2010;3:45. Jesmin J., Rashid M.S., Gene regulatory network reveals oxidative stress as the underlying molecular mechanism of type 2 diabetes and hypertension, BMC Med Genomics 3 (2010), 45.
 (12) Shamir R. Karlebach, G., Modelling and analysis of gene regulatory networks, Nature Reviews Molecular Cell Biology 9 (2008), 770–780.
 (13) R. Keller, Physical biology returns to morphogenesis, Science 338 (2012), no. 6104, 201–203.
 (14) L. Mendoza and E. AlvarezBuylla, Dynamics of the genetic regulatory network for arabidopsis thaliana flower morphogenesis., Theoretical Biology 193 (1998), no. 2, 307–319.
 (15) C. H. Waddington, The epigenotype, Endeavour 1 (1942), 18–20.