Forest structure in epigenetic landscapes

TL;DR

This paper introduces Epigenetic Forests, a novel tool for analyzing genetic regulatory networks to understand morphogenesis, demonstrated through modeling Arabidopsis flower development and accurately recovering its architecture.

Contribution

The paper presents Epigenetic Forests as a new method to study GRN structure and dynamics, providing insights into morphogenesis and cell fate determination.

Findings

Successfully modeled Arabidopsis flower development

Recovered flower architecture using genetic algorithm

Provided new insights into GRN structure during morphogenesis

Abstract

Morphogenesis is the biological process that causes the emergence and changes of patterns (tissues and organs) in living organisms. It is a robust, self-organising 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.

Tables

Table 1. Table 1: Equilibrium points: each component of the string q i subscript 𝑞 𝑖 q_{i} represents the state of a given gene (node): 0 is it is inactive and 1 if it is active. There are 13 characters in each q i subscript 𝑞 𝑖 q_{i} , one for each node.

Floral organ         Atractor
Inflorescence 1	$p_{I_1}=[0,0,0,1,0,0,0,0,1,0,0,0,0]$
Inflorescence 2	$p_{I_2}=[0,0,0,1,0,0,0,0,1,0,0,0,1]$
Inflorescence 3	$p_{I_3}=[0,0,0,1,0,0,1,0,1,0,0,0,0]$
Inflorescence 4	$p_{I_4}=[0,0,0,1,0,0,1,0,1,0,0,0,1]$
Sepals	$p_S=[0,1,1,0,1,1,0,0,0,0,1,0,0]$
Petals	$p_P=[0,1,1,0,1,1,0,0,0,1,1,1,1]$
Stamens	$p_T=[1,1,0,0,1,1,0,1,0,1,1,1,1]$
Carpels	$p_C=[1,1,0,0,1,1,0,1,0,1,1,0,0]$

Table 2. Table 2: Minimum energy chain: first column shows the state of each cell type in the chain, the second one the tree it belongs to and the third one its specialization energy.

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

Equations

$$c_S=152,c_P=160,c_T=3744,c_C=3608,$$

$$c_{I_1}=128,c_{I_2}=128,c_{I_3}=64,c_{I_4}=64.$$

$$C=[c_1,c_2,\dots c_L],$$

$$E_{SP}(c_1)=\sum _{i=1}^{n}d(c_i,c_{i+1})$$

$$E(C)=\sum _{i=1}^L{E}_{SP}(c_i).$$

$$I_i=[p_1,p_2,\dots ,p_{12}]$$

$$C_i=[c_S,c_1,\dots ,c_{13}]$$

$$fit(I_i)=\frac{\alpha }{E(C_i)}$$

$$I_i=[p_1,p_2,\dots ,p_{12}]$$

Full text

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.

††journal: Journal of Theoretical Biology

1 Introduction

Genetic Regulatory Networks (GRN) govern morphogenesis, the process responsible for producing the complex shape of full grown organisms from a few cells ([13], [8]) 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 ([5]).

In recent years GRN have gained a lot of attention and great advances have been made ([12], [6]). From cancer identification ([10]),([3]) and diabetes ([11]) to the formation of animals’ bodies ([4]) 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 ([2]), 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 ([9]), ([1]).

In ([2]), 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 ([15]), 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 ([14]) 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 $13$ nodes, each one corresponding to a specific gene 111The 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) [1]..

Each node has two possible states ([math] or $1$) 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 $N_{i}$ to node $N_{j}$ if node $Ni$’s state affects $N_{j}$ and a purple edge between node $N_{i}$ and $N_{j}$ if both nodes affect each other. The rules that the automata follows (obtained experimentally) can be found in ([14]).

There are $2^{13}$ possible initial conditions (binary numbers form [math] to $2^{13}-1$). The system is iterated, starting from each initial condition. It converges to eight different attractors 222For 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) [7]. Each equilibrium point has thirteen components (see table(1)).

The number of nodes (initial conditions land in each fixed point) in each tree will be the following:

[TABLE]

where the subscripts S, P, E and C correspond to sepals, petals, stamens and carpels respectively and $I_{i}$ to the inflorescences. ([7],[2]).

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 $S_{i}$ 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 $2^{13}$ nodes in the trees). There will be a directed edge from node $S_{i}$ to $S_{j}$ of the tree if (and only if), following the rules of the automata, the state of node $S_{i}$ leads to the one corresponding to node $S_{j}$ (see fig. (3)).

Proceeding this way until every state is analyzed, we obtain an in-tree (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 $8$ fixed states, we will end up with $8$ 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 $C$ of cell types. A chain will be formed by $L$ 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 $13$ entries (one for each gene) that can be in one of the two possible states we are considering ([math] or $1$). 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 $L$ links,

[TABLE]

where each link is a possible state of the automata, that is, $c_{i}=[x_{1},x_{2},\ldots,x_{13}]$ and $x_{i}=\{0,1\}$. The initial one, $c_{1}$ 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** (Hamming-distance)**

Given two strings of the same length $s_{i},s_{j}$, the Hamming distance between them $d(s_{i},s_{j})$ 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 $c_{i}$ to the root of the tree.

Definition 3.3** (Specialization path)**

The path from a given node $c_{i}$ in a tree to its root is called the specialization path of $c_{i}$ and we denote it as $SP(c_{i})$.

A chain specializes (differentiates) $C=[c_{1},c_{2},\ldots,c_{L}]$ when each of its links reaches a steady state (following the rules of the automata), that is if $c_{i}$ belongs to tree $T_{j}$, it must traverse its specialization path $SP(c_{i})$ reaching this way, the root of $T_{j}$.

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 $c_{1}$, $SP(c_{1})=\{c_{1}\to c_{2}\to\ldots\to c_{n}\}$ (where $c_{n}$ is the root of the tree we are on) is given by

[TABLE]

Definition 3.5** (Energy of a chain)**

The energy of a chain $C=[c_{1},c_{2},\ldots c_{L}]$ is given by

[TABLE]

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 $13$ links ($13$ being the number of components in each link) is redundant, so it will not be optimal (one with more than $13$ links can always be reduced to one with $13$). We can thus, limit our search to chains with $L=13$.

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 $c_{S}$ (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 $1$ it will become a [math] and viceversa). We can represent each individual as

[TABLE]

where $p_{i}\in[1,13]$ represents the component that will be modified in the following link, so that the chain $C_{i}$ generated by $I_{i}$ will be

[TABLE]

where $c_{1}(j)=c_{S}(j),j\neq p_{1}$, $c_{1}(p_{1})=mod(c_{S}(p1)+1,2)$ and in general $c_{i}(j)=c_{i-1}(j),j\neq p_{i}$ and $c_{i}(p_{i})=mod(c_{i-1}(p_{i})+1,2)$ (every component but one will be the same between neighbouring links, so that $d(c_{i},c_{i+1})=1$.

We start by randomly generating a population of $N$ individuals. $N$ 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 $C_{i}$, that is, we will compute its energy $E(C_{i})$ 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 $1$:

[TABLE]

where $\alpha$ 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 $\epsilon$.

4.2.1 Selection

Once we have a set of individuals (a population) $P=[I_{1},\ldots,I_{N}]$ 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 $1$. A random number $m\in[0,1]$ is chosen. The selected individual is the first one whose accumulated normalized value is greater than $m$. We repeat this process until $p$ 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 $I_{i}$ and $I_{j}$ using the selection operator. We will randomly choose two different crossover points $cp_{1},cp_{2}\in[1,L-1]$. The two individuals will be combined forming two new ones: $I_{i}^{n}$ and $I_{j}^{n}$. The first part of $I_{i}^{n}$, up to the crossover point $cp_{1}$ and from $cp_{2}$ to the end, will be copied from $I_{i}$, the second part (from $cp_{1}$ until $cp_{2}$), will be copied from $I_{j}$. Conversely, the first part of $I_{j}^{n}$ (up to $cp_{1}$) and from $cp_{2}$ to the end will be copied form $I_{j}$ and the second part from $I_{i}$ (see figure (14)). The two new individuals $I_{i}^{n}$ and $I_{j}^{n}$ each define new chains $C_{i}^{n}$ and $C_{j}^{n}$ 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 $80\%$.

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 $I_{i}$ using the selection operator. We will randomly choose two mutation points $mp_{1},mp_{2}\in[1,L-1]$. We will replace the value of the states $mp_{j}$ in $I_{i}$ by a randomly chosen one between $1$ and $13$ (see figure (15)). This will produce a new individual (with two mutated components):

[TABLE]

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 $20\%$ 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 $L=13$

• 2)

Maximum number of possible iterations allowed: $MaxG=200$

• 3)

Crossover probability: $cp=0.8$

• 4)

Mutation probability: $cp=0.2$

• 5)

Number of individuals per generation: $I=100$

• 6)

Elitism: $E=2$

After running the genetic algorithm an optimum solution is always found in less than $140$ iterations (generations). We show an execution where the minimum was found in $119$ iterations and the starting point is $0101100000000$ (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).

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.