Floral morphogenesis: stochastic explorations of a gene network epigenetic landscape

Abstract

In contrast to the classical view of development as a preprogrammed and deterministic process, recent studies have demonstrated that stochastic perturbations of highly non-linear systems may underlie the emergence and stability of biological patterns. Herein, we address the question of whether noise contributes to the generation of the stereotypical temporal pattern in gene expression during flower development. We modeled the regulatory network of organ identity genes in the Arabidopsis thaliana flower as a stochastic system. This network has previously been shown to converge to ten fixed-point attractors, each with gene expression arrays that characterize inflorescence cells and primordial cells of sepals, petals, stamens, and carpels. The network used is binary, and the logical rules that govern its dynamics are grounded in experimental evidence. We introduced different levels of uncertainty in the updating rules of the network. Interestingly, for a level of noise of around 0.5-10%, the system exhibited a sequence of transitions among attractors that mimics the sequence of gene activation configurations observed in real flowers. We also implemented the gene regulatory network as a continuous system using the Glass model of differential equations, that can be considered as a first approximation of kinetic-reaction equations, but which are not necessarily equivalent to the Boolean model. Interestingly, the Glass dynamics recover a temporal sequence of attractors, that is qualitatively similar, although not identical, to that obtained using the Boolean model. Thus, time ordering in the emergence of cell-fate patterns is not an artifact of synchronous updating in the Boolean model. Therefore, our model provides a novel explanation for the emergence and robustness of the ubiquitous temporal pattern of floral organ specification. It also constitutes a new approach to understanding morphogenesis, providing predictions on the population dynamics of cells with different genetic configurations during development.

Figure 1. Flower development and gene network underlying primordial floral organ cell-fate determination in Arabidopsis thaliana.

(A) The inflorescence meristem (IM in the Scanning Electron Micrography) is found at the apex of a reproductively mature plant. Within the IM, four regions can be distinguished. Interestingly, the experimentally observed gene activation configurations of each one of these regions are mimicked by the I1, I2, I3, and I4 attractors of the 15-gene GRN. Flower meristems arise in a helicoidal pattern from the flanks of the IM. The order in which floral meristems appear is indicated with numbers (1, oldest; 5, youngest). (B) Young flower meristems can be subdivided into four regions, each one containing the primordial cells that will eventually develop into the flower organs. In each floral meristem, the outermost region, which is first determined, will give rise to the sepal (se) primordium, the next to petals (pe) and finally, the primordial corresponding to stamens (st) and carpels (car) are determined in the center third and fourth whorls of the flower bud, respectively. (C) The mature flower of Arabidopsis thaliana. (D) I1, I2, I3, and I4 regions of the IM correspond to four of the attractors of the 15-gene GRN model. The expressed genes for each attractor are represented as gray circles, while the non-expressed genes correspond to white circles. (E) The other six attractors of the GRN model match gene expression profiles characteristic of sepal, petal (p1 and p2), stamen (st1 and st2), and carpel primordial cells. Black circles represent a gene (UFO) that can be either expressed or not expressed in the petal and stamen attractors, thus yielding two attractors for petal and stamen primordial cell-type. The gene activation profiles of the attractors recovered for the 15-gene GRN are congruent with the combinatorial activities of A, B, and C-type genes predicted by the ABC model of floral organ determination. See the Results section and , for details. (F) Gene regulatory network model underlying cell fate determination in the IM and the flower meristem. A-genes (red), B-genes (yellow), and C-genes (blue) from the ABC model are indicated in the network.

Figure 2. Heat map of the similarity matrix among the ten attractors of the GRN.

A strict consensus phenogram was obtained for the GRN attractors (vectors of zeros and ones) by using the Manhattan distance similarity index (see Methods). This phenogram is shown below the attractors that are ordered along the X and Y axes of the heat map. Attractors that group together had the highest similarity indexes between them (i.e. the lowest Manhattan distance). Color scale: darker colors indicate more similar, while lighter ones indicate more different attractors in the pairs compared.

Figure 3. Temporal sequence of cell-fate attainment patterns under the Boolean dynamics with noise.

Maximum relative probability (“Y” axis) of attaining each attractor, as a function of iteration number or time (“X” axis). (A) Probability of attaining each attractor (i.e., cell type) obtained by multiplying the Markov matrix M by a population vector initialized at the sepal attractor. The error probability in computing this graph was η = 0.03. The most probable sequence of cell attainment is: Sepals, petals, carpels, and stamens. (B) Probability of attaining each attractor (i.e., cell type) at each iteration when 80000 randomly chosen “sepal” configurations were selected and followed for 140 steps. Noise was introduced in the updating of each gene independently, with a η = 0.03 probability at each iteration. The probabilities for the petal (p) and stamen (st) attractors correspond to the sum of p1+p2 and st1+st2, respectively. All maxima correspond to 100 because each absolute probability value was divided by the maximum of each attractor's curve (see Results and Methods). Equivalent graphs to those in (A) and (B) for η = 0.01 are shown in (C) and (D), respectively.

Figure 4. Changes in the basins of attraction of the continuous model with respect to the Boolean model.

Color map of the probability P(n|m) that a microscopic configuration whose associated Boolean configuration belongs to the basin of attraction of attractor m, ends up in attractor n using Glass dynamics. Note that the main transitions occur along the diagonal where attractors are reached by both dynamics (Boolean and Glass); however, the non-diagonal elements indicate that two microscopic configurations that correspond to the same Boolean configuration may end up in different attractors.

Figure 5. Effects of the choice of the relaxation time on Glass dynamics with noise.

Two typical realizations of Glass dynamics for a given gene x_n showing that the choices of the relaxation time τ and the perturbation time Δt_p do not affect the qualitative dynamics, so long as Δt_p>τ. Both trajectories started from the same initial conditions, and were followed through the same set of perturbations. The black trajectory corresponds to Δt_p = 2.5 and τ = 1, whereas the red trajectory corresponds to Δt_p = 2.5 and τ = 1/20.

Figure 6. Temporal sequence of cell-fate attainment patterns under the Glass dynamics with noise.

Maximum relative probability (“Y” axis) of attaining each attractor as a function of iteration number or time (“X” axis). (A) The maxima of the cell-fate curves are attained in a particular sequence in time, which in this case is sepal, petal, stamen, and carpel. Parameters used: dt = 0.01, τ = 1, and Δt_p = 2.5. (B) When the simulations mimic the Boolean case (dt = 1, τ = 1 and Δt_p = 1; see Results and Methods), a temporal pattern identical to that of the Boolean dynamics was obtained, with a sequence of sepal, petal, carpel and stamen. The noise used in both cases was η = 0.03. Although the Boolean and Glass dynamics need not coincide in general, for the case of the A. thaliana GRN, both models provide similar predictions. Simulations show that the order of emergence of the stamen and carpel maxima, as compared to the Boolean model, may depend on the precise values of the kinetic constants.

Figure 7. Schematic representation of the epigenetic landscape generated by a stochastic exploration of the GRN for flower development.

This schematic landscape is equivalent to the Epigenetic Landscape proposed by C.H. Waddington (1957). Basins comprise the cell genetic configurations that lead to attractors (in this case, gene arrays characteristic of floral organ primordial cell-types: Sepals, petals, stamens, and carpels. See Figure 1 and Discussion). Each cell fate is associated to the GRN configuration corresponding to each of the attractors. The arrows represent transitions among attractors. The transition from inflorescence to sepal attractor might be biased or determined by an inducer. The numbers associated to the arrows represent the sequence of transitions among attractors: From sepals to petals, and then to carpels and stamens.