Abstract
Background
The genetic control of floral organ specification is currently being investigated by various approaches, both experimentally and through modeling. Models and simulations have mostly involved boolean or related methods, and so far a quantitative, continuoustime approach has not been explored.
Results
We propose an ordinary differential equation (ODE) model that describes the gene expression dynamics of a gene regulatory network that controls floral organ formation in the model plant Arabidopsis thaliana. In this model, the dimerization of MADSbox transcription factors is incorporated explicitly. The unknown parameters are estimated from (known) experimental expression data. The model is validated by simulation studies of known mutant plants.
Conclusions
The proposed model gives realistic predictions with respect to independent mutation data. A simulation study is carried out to predict the effects of a new type of mutation that has so far not been made in Arabidopsis, but that could be used as a severe test of the validity of the model. According to our predictions, the role of dimers is surprisingly important. Moreover, the functional loss of any dimer leads to one or more phenotypic alterations.
Background
For various plant species, floral organ specification has been successfully linked to spatial gene expression patterns according to the wellknown ABC model [1]. This model has recently been extended to five gene classes (ABCDE) to explain novel floral mutants and to accommodate functions that specify ovule development and the establishment of floral context [27]. Despite these modifications, however, the ABCDEmodel remains a static, qualitative model that does not describe the detailed molecular interactions involved, nor the temporal and spatial gene expression patterns that these interactions induce.
To model the molecular interactions involved in floral organ formation, various approaches have been used, mostly in terms of boolean networks. A boolean network approach was successfully applied to recover known stable states and to predict the existence of unknown interactions, [810]. Also, a stochastic type of boolean network, and a differential equations model, that can be considered as a first approximation of kineticreaction equations, have been proposed [11]. These types of models are especially suited for qualitative analysis of large model structures. The validity of a candidate model can be tested by comparing the steady states of the model with those measured experimentally. In [12], a general review is given on the various modeling approaches applied to gene regulatory networks, ranging from basic logical models to very extensive stochastic simulation algorithms, and a review specifically on stochastic methods is given in [13]. In [14], a stochastic model of the autoregulatory loop of the Btype genes PISTILLATA (PI) and APETALA (AP3) in Arabidopsis and Antirrhinum is described.
Ordinary differential equations (ODE) have been used extensively to model gene regulatory networks, including the Notch signaling pathway [15], the Zebrafish Somitogenesis Clock [16], the carbon starvation response in E. coli [17], and the toggleswitch gene network [18]. This type of model allows a quantitative, continuoustime analysis. However, for quantitative reliability, detailed parameter information is essential. This information is often not available, and instead the parameters are estimated indirectly by an identification procedure. In [19] an overview is presented of ODE networks that have been identified or are suitable as benchmark test.
Recently, gene expression data sets have become available for the genes involved in specification of floral organ identity. In [20,21], time series of gene expression are presented for each class of genes in the ABCDE group. For most ABCDE genes, the majority of which are members of the MADSbox transcription factor family, the way in which they are activated is known from experiments [22,23]. Furthermore, it has been shown that MADS box transcription factors regulate their own and each other's expression via various autoregulatory loops (reviewed in [24,25]). These two information sources open the door for ODE model development. There is also considerable evidence that MADS proteins play a crucial role in the autoregulatory repression or activation of specific sets of target genes [2528].
The consensus MADS domain protein target site (CArGbox) is a palindromic sequence and structural analysis of the SRF MADS domain region bound to DNA has revealed its binding in a dimeric form [29]. Based on the identification of higherorder complexes and the fact that more than two different MADS domain proteins are essential for specifying organ identity, it has been hypothesized that plant MADS transcriptional complexes consist of quartets assembled from two active dimers [3033]. The structure of this paper is as follows. First, we develop an ODE model for MADS box gene expression. The role of dimers is explicitly incorporated in this model. Second, we show that this model is able to simulate the dynamics of experimental time series data. Third, after the model has been fitted to experimental data, the predictive power of the model is assessed by comparing model predictions of ABCDE gene mutants with independent mutant experiments. Finally, we use the model to predict the effects of changes in the topology of the underlying proteininteraction network. We conclude that the model has a good predictive power with respect to mutations. A simulation study of a new mutant type in which the formation of specific dimers is disrupted, shows that each dimer function is essential for proper organ formation.
The ABC model for Arabidopsis
The ABC model links spatial gene expression patterns to phenotypes. Figure 1 shows the expression domains of the five gene classes corresponding to the four types of floral organs: sepals, petals, stamens and carpels (including ovules). For clarity, the carpels and ovules will here be regarded as one organ type. The figure illustrates that, for example, sepal identity is determined by (high) expression of A and E type genes, and this normally occurs in whorl 1, the most outer whorl of the floral meristem. In Arabidopsis, the five gene classes in the ABCDE model comprise several redundant genes. The A type genes are represented by APETALA1 (AP1) and APETALA2 (AP2), the B type by AP3 and PI, the C type by AGAMOUS (AG), the D type by SHATTERPROOF1 and SHATTERPROOF2 (SHP1, SHP2) and SEEDSTICK (STK), and the E type by SEPALLATA1SEPALLATA4 (SEP1, SEP2, SEP3 and SEP4). SEP13 are expressed in whorl 24, and SEP4 is expressed in all whorls. The A class gene AP2 is exceptional in that it is the only floral organ identity gene that does not belong to the MADS domain transcription factor family. Although many of the ABCDE genes have been studied extensively, the detailed genetic and molecular interactions among the various redundant genes are still not known comprehensively, and for a number of genes sufficiently detailed expression data are lacking. Therefore, as a first model simplification, we only take MADS box genes into account and we assume that (partially) redundant genes have similar dynamic expression patterns and similar interactions, and can therefore be represented by a single gene from each clade. This is a common approach in literature in dealing with redundancy [8,9,11,14]. The following representative genes are selected to represent the five ABCDE functions: AP1 (A), AP3 (B), PI (B), AG (C), SHP1 (D) and SEP (E). Here, SEP denotes SEP3 in whorl 24, and SEP1SEP4 in whorl 14.
Figure 1. The ABCDE model for flower organ determination in Arabidopsis. The figure has to be read columnwise. E.g., in the sepalwhorl 1, genes A and E are dominantly expressed.
Methods
Network properties
Dimers are known to play an important role in the dynamics of the MADS protein network and represent the minimal structural unit essential for DNA binding [27,29]. We therefore explicitly include dimers in our model and consider their regulatory functions as known or suggested by indirect genetic evidence (Table 1) [25,34,26]. Figure 2 shows the network interactions graphically. Table 2 gives the expected protein expression patterns in the different whorls from day 2 to 5 of meristem development. Before day 2, no initiation of floral organ primordia and differentiation of the different floral organs takes place, and the genes are expressed uniformly over the floral meristem.
Table 1. Dimers that are supposed to be repressing (right column) or activating (middle column) the transcription of the individual MADS domain proteins (left column).
Table 2. Expressed MADS box genes for each floral whorl from day 2 to 5 of meristem development. Day 2 means the end of the second day. SEP (E) is abundant everywhere at all times.
Figure 2. Graphical representation of the interactions in Table 1.
Dynamical model
To model the dynamical behavior of the MADS box genes, we write down the governing differential equations, which are based on the following set of assumptions: (a) Intercellular diffusion of MADS box proteins is, on average, small. With diffusion ignored, ODE's instead of partial differential equations (PDE's) can be used. The different whorls are modeled with whorlspecific triggering mechanisms, representing the timed initial activation of the MADS domain proteins by nonMADS factors, such as LEAFY (LFY), and WUSCHEL (WUS) [22,23]. (b) After activation, MADS box gene transcription is regulated by autoregulatory mechanisms in which protein dimers play an essential role. (c) Transcription only occurs when one or more activation sites on the DNA are occupied, and simultaneously all repression sites are empty. This assumption is commonly made in modeling of gene regulatory networks [35]. (d) The delay effect of translation is neglected, i.e. transcription immediately leads to protein formation. (e) Dimer decay into nonfunctional components is small compared to decay into functional monomers [36,37] and is therefore ignored. (f) During the first five days of meristem development, the average cell size remains approximately the same [38].
The dynamics of the dimer concentrations consists of the association rate of monomers into dimers, minus the dissociation rate of dimers into monomers. Denoting by x_{i }the concentration of monomer i and by [x_{i}x_{j}] the concentration of the dimer of proteins i and j, we have the following balance equation
The proteins represented by the variable x_{i }are listed in Table 3.
Table 3. Notation for the variables that represent the MADS proteins.
The dynamics of monomer concentrations is more complex. It depends not only on the dimer association/dissociation rates, but also on transcriptional regulation and decay into nonfunctional elements. Transcriptional regulation is modeled by MichaelisMenten functions, in which β represents the maximum transcription rate, Km the halfmaximal activation or repression, and dc the decay rate. As additional elements to the model we include two triggers, p_{2}(t, w) and p_{4}(t, w), which govern the expression of genes AP3 and AG, respectively, at time t in whorl w. The expression of AP3 is regulated by the genes LFY and UNUSUAL FLORAL ORGANS (UFO) [39], and the expression of AG is regulated by the genes WUS and LFY [40,41]. Since these terms are the only timeand whorldependent components in the model, they are responsible for cell differentiation. The triggers are essential to drive the network into four different steady states, where each one corresponds to a different organ identity. The biological mechanism that is responsible for the trigger is not modeled here. Quantitative information on the amount of protein generated by the triggers is not available, but their timing is known. Because autoregulatory loops can maintain the expression of the MADS box genes after induction, the duration of the triggers is set to one day. The triggers p are thus assumed to take on a nonzero constant value between day 1 and 2 only, and otherwise are set to zero. Altogether, this gives the following model for the monomer dynamics:
Here, the first fractions on the right hand sides denote activation or repression by MichaelisMenten kinetics, followed by a decay term. The last terms denote the rates of dimerization, which, when positive, act to decrease monomer concentrations. Per whorl, the network dynamics is governed by equation sets (1)(2), which involves 13 equations, 13 variables, and 51 parameters. To enable parameter estimation, we reformulate the model entirely in terms of monomer concentrations. Another advantage to this is that elimination of the dimer variables and dynamics considerably simplifies the analysis, since it reduces the number of equations and variables to 6. This makes the search algorithm easier to implement and faster. Reformulation is done, first, by applying a time scale decomposition to (1). For a comprehensive background to this technique, see e.g. [42], p.168. The time constant is days, which implies that on a scale of days, the dynamics of dimer formation is very fast. This justifies the use of a quasisteady state approach, in which the dimer concentrations are fully determined by the instantaneous monomer concentrations. This effectively comes down to setting the time derivatives in (1) to zero. By doing so, the dimer equations (1) take the form
and these are inserted into (2), using the chain rule:
If these expressions are substituted in (2), we obtain a system of 6 equations with 6 variables:
Here, The ordering of γ_{j }is discussed below. We are aware that these equations now attain a form that is quite unusual for ODE's, since the derivatives are also present in the right hand sides. We explain in the following section why this form is still useful in the context of parameter estimation.
Parameter estimation
The 37 parameters in (5) need to be estimated from measured time series. This estimation is done in three steps. First, the allowed parameter ranges are defined. Second, a data set is presented and converted into the desired whorlspecific form. Third, the network and data set are decoupled to allow a successful identification procedure.
Parameter values
Because the number of parameters is relatively large compared to the number of data points, straightforward estimation might be problematic. Hence, we choose to treat different parameters on a somewhat different footing, depending on available biological knowledge about allowed parameter ranges.
In [43] values are given for β in the range of [3, 41] nMmin^{1 }(α, Table 3). We therefore take [1, 50] nMmin^{1 }as a reasonable range for β. In the same table in [43], Km (there named THR) is given in the range of [10^{2}, 10^{3}] nM, and in [44] values used for Km (there named )) are in the range of [10, 10^{2}] nM. Therefore, a reasonable range for β is [1, 50] nMmin^{1}. A range for decay of dc ∈ [10^{3}, 10^{1}] min^{1 }is given in [45]. Ranges for association and dissociation constants are K_{on }∈ [10^{3}, 1] nM^{1}min^{1 }and K_{off }∈ [10^{3}, 10] min^{1}[45]. The relative interaction strengths between dimerising proteins are based on expert knowledge:
where {x_{i} x_{j}} denotes the value of K_{on}/K_{off }corresponding to the dimer [x_{i}x_{j}]. Based on this information, we fix the values of these parameters at K_{off }= 1, γ_{1 }= 1, γ_{2 }= 10^{1}, γ_{3 }= 10^{1}, γ_{4 }= 10^{2}, γ_{5 }= 10^{2}, γ_{6 }= 5.10^{3}, and γ_{7 }= 5.10^{3}nM^{1}.
Data manipulation
Data from [21] contain the mRNA signal intensities of the six genes included in our model, at the first five consecutive days of floral meristem development. AP1, which is normally activated by LFY and FLOWERING LOCUS T (FT) [39,46], is induced here artificially at time point 0. The measured SEP3 concentrations are expected to be representative for the concentrations of SEP1SEP4 in whorl 24, and for SEP4 in whorl 1. The intensities are averages over the whole meristem. Since we need whorlspecific data, the data set is transformed from average intensities to whorl specific protein concentrations in five steps.
1. The data set is scaled uniformly from mRNA intensity to protein concentration, such that the average concentration is 10^{3}nM (in [45] a range of 10^{2}10^{4 }nM is mentioned for transcription factors in eukaryotic cells). Here we implicitly assume that the microarray signal intensities have a linear correspondence to the protein concentrations. This is based on the observation that spatial mRNA expression levels and protein levels correspond well to each other for some of the ABC class MADS transcription factors [47].
2. If in a whorl a gene is not expressed, the protein concentration is set to 1% of the value of a protein that is expressed. This is based on visual interpretation of confocal images from [47].
3. The gene expressions per time point and whorl are based on [20] and are given in Table 2.
4. The relative whorl sizes are obtained from confocal images from [47]. From day 2 to day 5, organ identity comes into play. The shapes and relative sizes stay approximately the same between day 2 and 5. At the end of day 2, the sepals have a volume of 1.1·10^{4}μm^{3}, the petals of 2.7·10^{4}μm^{3}, the stamens of 2.9·10^{4}μm^{3}, and the carpels of 1.1·10^{4}μm^{3}.
5. The mass balance for the concentration of protein i is
Here, w runs over the whorls, is the average concentration of protein i from the data set, the concentration of protein i in organ w, and V^{w}(t) the organ volume. Further,
with the percentage of expression (1% when there is no expression, 100% when the gene is expressed), and α_{i}(t) a scaling factor. To determine α_{i}(t), equation (8) is inserted into (7), which yields the expression
Network decoupling
With the γ_{i }values given, the system of ODE's (5) still contains 37 parameters that need to be estimated. This puts a high computational demand on the search algorithm, which we propose to alleviate by using a decoupling procedure [4851]. This approach is based on a simple, highly effective idea. Let us explain this for the decoupling of equation 5(a), which has the form:
with par the set of parameters in this equation. For concentrations x_{2, }.., x_{6 }and we take the data and interpolate them with a forward Euler method. This basic interpolation scheme is straightforward and does not introduce substantial interpolation errors. In the end the decoupled network of monomers has the same quality of data fit as the coupled dimer network, Figures 3 and 4.
Figure 3. Simulated dynamics of the decoupled model (5) (solid lines) of the monomers, together with the data points for the four organs.
Figure 4. Simulated dynamics of the coupled network (1)(2) (solid lines) of the concentrations of proteins that are part of dimer complexes, together with the data points for the four organs.
The resulting functions x_{2}(t), ..., x_{6}(t) and are substituted in (10) so that in the resulting equation x_{1 }is the only variable. This equation is integrated and by fitting the calculated values of x_{1 }to the data for x_{1}, the parameters par are found. This procedure thus leads to estimates for the subset of parameters in (10). Similarly, the other parameters in 5(b)5(f) are estimated by decoupling the equation under consideration from the others. Note that this reduction method is applicable thanks to the fact that no parameter appears in more than one equation.
The measured concentrations are those of the total amount of x_{j }, both in monomer and dimer form, which are denoted by . To calculate the monomer concentrations from the data, we use mass balances. From the dimers listed in section Network properties, we find that
The quasisteady state equations (3) are inserted into (11), which yields a set of nonlinear equations in the monomer variables. Finally, the monomer concentrations are estimated by a nonlinear search algorithm for each time point and whorl in the data set.
Parameter identification
The identification criterion is to minimize the least squares error
by optimization over the parameter vector a, that consists of the unknown parameters in (5). Here are the data points in whorl w for protein j at time t_{i }(in days), and x_{j, w}(t_{i}, a) the concentrations that are predicted by our model for some choice of parameter vector a. The optimization of a is carried out by the Matlab routine "lsqnonlin", which is a gradientbased search method [52,53]. The initial concentrations (at i = 0) are taken equal to the data points. For the integration we used the Matlab ode23 DormandPrince algorithm [54].
As for the choice of initial values, we investigated several strategies and found that the following choice led to the fastest convergence of the search algorithm. The initial parameter values K_{m }were chosen such that for typical concentration values the MichaelisMenten functions attain their maximal slopes and therefore are highly sensitive for parameter variations.
The search space for a is confined as much as possible to the parameter ranges that are listed in section Parameter values. In Additional file 1 a robustness analysis is presented which aimed to assess whether the optimal parameter values that are retrieved are sensitive to the choice of the γ's. It turned out that the local minima are robust against variation in any γ by at least a factor two.
Additional file 1. Sensitivity analysis of the optimal parameters with respect to binding affinities.
Format: PDF Size: 93KB Download file
This file can be viewed with: Adobe Acrobat Reader
Results
Parameter estimation
The estimated parameter values are listed in Table 4.
Table 4. Identified parameters for model (5). β_{i,j }is the maximal transcription rate in nM day^{1 }for the j th MichaelisMenten function of gene i, Km_{i,j }the corresponding halfmaximal activation in nM, d_{i }the decay in day^{1}, and P_{i }the trigger in nM day^{1}.
According to section Parameter values, the maximal transcription rate, which is the sum of the β's in each equation, should lie in the range of [1.4·10^{3}, 7·10^{4}] nM/day. All the (sums of the) β's are within this range, except for β_{5,1}, which is a factor 3.5 lower. The values of K are all within the range of [10, 10^{3}] nM^{1}, with the exception of Km_{2,2 }and Km_{4,3}, which are only slightly higher. All decay rates are within the range of [1, 1.5·10^{2}] day^{1}, except dc_{4}, which is somewhat higher.
Figure 3 shows the simulated decoupled dynamics of the monomer proteins in the identified model (5), together with the data points.
For convenience, highly expressed genes will be called "on" and lowly expressed genes will be called "off". AP1 has a good fit in the sepals and petals, where the gene is on. In stamens and carpels, AP1 is a factor 10 too high from day 25, but still a factor 5 lower than the onlevel. The contributions of the triggers of AP3 (petals and stamens) and AG (stamens and carpels) between day 1 and 2 are clearly distinguishable. PI shows an overshoot between day 1 and 2 in the sepals and carpels. It has to be mentioned that the logscale visually magnifies errors in the low concentration regime. Since the assumption that the offgenes have 1% of the concentration of the ongenes is an estimation that induces unavoidable errors in the low regimes, these are considered less important. A parameter set is found that fits the monomeric concentration data quite reasonably. It is well known [30] that cell identity is determined by complexes of the ABCDE genes, instead of only by the monomers. Therefore, the monomer concentrations of the model (1)(2) are converted into dimer concentrations (i.e. concentrations of proteins that are part of a dimer), using (3) and the mass balance (11). The data set is converted similarly. Figure 4 shows the profiles of the concentrations of proteins that are part of dimer complexes in the resulting coupled network. Since AP3 and PI only form the dimer [AP3 PI], their concentrations are equal. Therefore AP3 is omitted in the dimer concentration plots from now on. The simulated concentrations remain more or less consistent after day 5, which one would expect. It could be expected that the fit of the coupled network in Figure 4 is less accurate than the fit of the decoupled equations in Figure 3. Conversion into dimer concentrations could result in further inaccuracies, since the parameters are fitted on monomeric data. However, comparison with Figure 3 shows that there are not many substantial differences between the coupled and decoupled network. Figure 5 shows for each whorl and each gene the mean relative error , which is defined as follows
Figure 5. Mean relative error between simulation and experiment as defined in equation (13). The horizontal axis corresponds to the variables [x_{1, }.., x_{6}].
where x is the simulated and x_{d }the measured gene expression. In (13) we take on purpose not the squares of the differences, since epsilon should measure whether x_{d }is systematically overor underestimated over five days. Moreover, the differences xx_{d }have been normalized by x_{d }in order to make a comparison between concentrations of different orders of magnitude possible. Based on Figures 4 and 5 we make the following observations:
From a qualitative point of view, the fitted model reproduces the expression data very reasonably. This implies that the topology of the model apparently includes the most relevant interactions. AP1 has a reasonable fit in the sepals and petals, where it is on. In the stamens and carpels it is too high after day 2, by a factor 8. However, the AP1 values are still a factor 12 below the onlevel. AG has a good fit for the onlevels in stamens and carpels, and in the sepals and petals it compromises between the high values at day 1, and the low levels later on. SHP is on only in the carpels after day 3, and here the fit is accurate. In the other organs there is a compromise between the high levels during day 03, and the low levels at days 45. SEP has an accurate fit for all organs.
Model validation by mutants
A powerful method to test the validity of our model is to compare the expression measurements in plants with mutations in particular MADS transcription factors with the MADS protein concentrations predicted by the model. The mutations can be simulated by either setting the initial concentration and production of protein from a knocked out gene equal to zero, or by fixing the concentration to a high level if the gene is ectopically expressed. We applied this test for five known mutants in Arabidopsis. For four mutants, the model predicted the correct phenotypes very well. One mutant was only partly predicted correctly. This confirms the predictive power with respect to genetic mutations.
In the first mutant, AP3 is missing. According to [55], the second whorl grows sepals, and the third whorl grows carpels. Figure 6 shows that our model predictions are in agreement with these phenotypic alterations. Indeed, we observe that the expressions in the second whorl agree with the first, and therefore they develop the same organs. The same applies to the third and fourth whorl. The expression levels of the Bgene PI become so small that they are not visible in Figure 6, due to the logarithmic scale. In the stamens the Bgenes are off, and since in this model SHP is repressed by the AP3PI dimer, the third whorl develops the same SHP levels as the fourth whorl. This repression assumption has never been proven, but is now supported by the outcomes of this simulation. The results of the remaining mutants are mentioned in brief hereafter. A more elaborate discussion of the simulation results is presented in Additional file 2.
Figure 6. Simulated dynamics for the AP3 = 0 mutant: the second whorls grow sepals, and the third whorl grows carpels. The data points denote the wildtype expression levels, and they are shown to compare the mutant dynamics with.
Additional file 2. Discussion on the simulation results of the mutants.
Format: PDF Size: 157KB Download file
This file can be viewed with: Adobe Acrobat Reader
In the second mutant, PI is missing. The same phenotype occurs as with AP3 missing [1], and this is in agreement with the model predictions. In the third mutant, AP3 is ectopically expressed. According to the model prediction, the fourth whorl organs have stamen expression, and the first whorl organs have petal identities. However, according to [36], the fourth whorl organs develop as stamens, but there is no change in the identities of the first whorl organs. The model prediction is thus only partially correct. In the fourth mutant, AG is ectopically expressed. According to the model prediction, the first whorl grows carpels, and the second whorl grows stamens. This is in agreement with the findings in [56]. In the fifth mutant, AG is absent in all whorls. According to [57], AG is involved in terminating the meristem activity in the flower after the formation of the fourth whorl carpels. Without AG, only sepals and petals are formed and the inner part of the meristem starts to develop new floral buds that reiterate the formation of sepals and petals ad infinitum. The AG stop mechanism is not modeled here, so in our simulations only 4 whorls of floral organs develop. Without the C gene, the A gene is not repressed in any floral organ, and will therefore be expressed in the whole floral meristem throughout flower development. According to the ABC model in Figure 1, whorls 14 will have sepalspetalspetalssepals identity, respectively, and this is exactly what our model predicts.
Model predictions
Recently we developed a bioinformatics method to predict specific sites in protein sequences where mutagenesis changes protein interaction specificity [58]. For various MADS proteins, interaction specificity was changed experimentally, and validated using yeasttwohybrid experiments. This is essentially a new type of mutation experiment, for which the phenotypic effects have not yet been investigated in planta. This kind of experiments can however easily be mimicked with the present model. We present results of this type of simulations to allow for the possibility of comparing them with future experiments. As a model prediction, the interaction strength of each dimer γ_{i }is set to zero, one at a time. This comes down to removing one specific dimer from the network. Table 5 shows the predicted phenotypic alterations.
Table 5. Predicted organ mutations for each removed dimer binding affinity.
It turns out that all dimer mutations show very clear organ conversions. It is interesting to see that all dimers play an important role in organ specification, even the dimers with very low association rates. Experimental verification of these predictions can provide a valuable tool for a more accurate determination of parameter values, as well as model structure, without the need for costly time series of expression data.
Discussion
The mutant simulation experiments show realistic results, and only one mutant (out of five) could not fully be reproduced. This indicates that we developed a very useful model, despite that it was based on only a limited amount of quantitative data. All seven dimer mutants are predicted to have phenotypic effects. In five out of seven mutations, double organ conversions occur, and for one mutation no floral organs formed at all, as expected in view of the known importance of the SEP proteins for determining floral organ identity.
These predictions, however, have to be interpreted with caution. The simulated dynamics depends on parameter values that have an uncertainty intrinsic to the sparse and uncertain data set used. Hence, accurate quantitative predictions are yet out of reach. Nevertheless, on a qualitative scale the mutant predictions suggest that the functional loss of each dimer leads to phenotypic alteration, and that therefore each dimer plays an essential role in the regulatory network. At this moment, our group is performing confocal image analysis of MADS protein expression patterns and levels to obtain a highquality and quantitative protein expression data set. We are also investigating the existence of additional genetic interactions. We expect that this will result in a complete network with accurate parameter values, that opens the door for testing other candidate model structures. For example, an additional hypothesis on the mode of action of the MADS proteins in the gene regulatory network that controls floral organ specification is proposed in [31] and states that the MADS proteins act in quaternary complexes, consisting of two dimers. Over the last few years more and more evidence has become available showing the formation [59,30,32] and DNA binding capacity of such complexes [59,60].
Conclusions
We proposed an ODE model for the dynamics of six genes that regulate floral organ identity in Arabidopsis. The model describes transcription regulation, mass balance, dimer formation, decay and cellfate determining trigger mechanisms. The parameters are estimated by an identification method that comprised a network decoupling method. The data set used consists of microarray intensities from the literature. The resulting model is validated by predicting the phenotypes of five mutants known from literature. Also, some new model predictions are made for an in vitro type of mutants, in which the formation of specific dimers is artificially repressed. Thanks to its welldefined mathematical and biological foundation, the model can be easily extended with additional biological knowledge. The model structure allows a decoupling procedure that seems to be a promising identification technique. Its application is apparently generic for ODE models of gene regulatory networks.
Experimental verification of the dimer mutants could provide a valuable tool for confirmation, falsification and/or model refinement. In systems biology, the cyclic interaction between 'wetlab' experiments and 'drylab' modeling plays a defining role. In this light we have presented here several mutants for which our model predicts a phenotype, in the hope to inspire biologists to test and study these special cases experimentally.
Authors' contributions
SM: Mathematical model analysis, simulation studies and software development, wrote the article. AD, RI, KK, GA and RH: Model development, data analysis and manipulation, scientific background. MG and JM: Mathematical model analysis and cowrote the article. All authors read and approved the manuscript.
Acknowledgements
We would like to thank Karel Keesman for useful discussions on model identification, and Susan Urbanus for providing us with the necessary confocal image data. This work was financed by the Netherlands Consortium for Systems Biology (NCSB) which is part of the Netherlands Genomics Initiative/Netherlands Organisation for Scientific Research and the Netherlands Organization for Scientific Research (NWO, NWOVENI Grant 863.08.027 to A. van Dijk).
References

Coen E, Meyerowitz E: The war of the whorls: genetic interactions controlling flower development.
Nature 1991, 353:3137. PubMed Abstract  Publisher Full Text

Ferrario S, Immink R, Angenent G: Conservation and diversity in flower land.
Current opinion in plant biology 2004, 7:8491. PubMed Abstract  Publisher Full Text

Causier B, SchwarzSommer Z, Davies B: Floral organ identity: 20 years of ABCs. [http:/ / www.sciencedirect.com/ science/ article/ B6WX04XKBYX02/ 2/ 72df2c458605237d76c5635658302624] webcite
Seminars in Cell & Developmental Biology 2010, 21:7379.
[TumorStroma Interactions; Flower Development].
PubMed Abstract  Publisher Full Text 
Colombo L, Franken J, Koetje E, vanWent J, Dons HJM, Angenent GC, van Tunen AJ: The Petunia MADS Box Gene FBP11 Determines Ovule Identity.
The Plant Cell 1995, 7(11):18591868. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Pelaz S, Ditta GS, Baumann E, Wisman E, Yanofsky MF: B and C floral organ identity functions require SEPALLATA MADSbox genes.
Nature 2000, 405:200203. PubMed Abstract  Publisher Full Text

Favaro R, Pinyopich A, Battaglia R, Kooiker M, Borghi L, Ditta G, Yanofsky MF, Kater MM, Colombo L: MADSbox protein complexes control carpel and ovule development in Arabidopsis.
Plant Cell 2003, 15:26032611. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Krizek BA, Fletcher JC: Molecular mechanisms of flower development: an armchair guide.
Nat Rev Genet 2005, 6:688698. PubMed Abstract  Publisher Full Text

EspinosaSoto C, PadillaLongoria P, Buylla E: A gene regulatory network model for cellfate determination during Arabidopsis thaliana flower development that is robust and recovers experimental gene expression profiles.
The Plant Cell 2004, 16:29232939. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mendoza L, Thieffry D, AlvarezBuylla ER: Genetic control of flower morphogenesis in Arabidopsis thaliana: a logical analysis.
Bioinformatics 1999, 15:593606. PubMed Abstract  Publisher Full Text

AlvarezBuylla ER, Azpeitia E, Barrio R, Bentez M, PadillaLongoria P: From ABC genes to regulatory networks, epigenetic landscapes and flower morphogenesis: Making biological sense of theoretical approaches. [http:/ / www.sciencedirect.com/ science/ article/ B6WX04XP381V6/ 2/ 5e0534192b4f9a12b9f6a54f2cb94735] webcite
Seminars in Cell & Developmental Biology 2010, 21:108117.
[TumorStroma Interactions; Flower Development].
PubMed Abstract  Publisher Full Text 
AlvarezBuylla ER, Chaos A, Aldana M, Bentez M, CortesPoza Y, EspinosaSoto C, Hartasnchez DA, Lotto RB, Malkin D, Escalera Santos GJ, PadillaLongoria P: Floral morphogenesis: stochastic explorations of a gene network epigenetic landscape.
PLoS ONE 2008, 3:e3626. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Karlebach G, Shamir R: Modelling and analysis of gene regulatory networks.

Pahle J: Biochemical simulations: stochastic, approximate stochastic and hybrid approaches.
Brief Bioinformatics 2009, 10:5364. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lenser T, Theissen G, Dittrich P: Developmental robustness by obligate interaction of class B floral homeotic genes and proteins.
PLoS Comput Biol 2009, 5:e1000264. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Agrawal S, Archer C, Schaffer DV: Computational models of the Notch network elucidate mechanisms of contextdependent signaling.
PLoS Comput Biol 2009, 5:e1000390. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Cinquin O: Repressor dimerization in the zebrafish somitogenesis clock.
PLoS Comput Biol 2007, 3:e32. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Porreca R, Drulhe S, de Jong H, FerrariTrecate G: Structural identification of piecewiselinear models of genetic regulatory networks.
J Comput Biol 2008, 15:13651380. PubMed Abstract  Publisher Full Text

Farcot E, Gouzé JL: A mathematical framework for the control of piecewiseaffine models of gene networks. [http://wwwsop.inria.fr/virtualplants/Publications/2008/FG08] webcite
Automatica 2008, 44(9):23262332. Publisher Full Text

Gennemark P, Wedelin D: Benchmarks for identification of ordinary differential equations from time series data.
Bioinformatics 2009, 25(6):780786. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Schmid M, Davison TS, Henz SR, Pape UJ, Demar M, Vingron M, Schlkopf B, Weigel D, Lohmann JU: A gene expression map of Arabidopsis thaliana development.
Nat Genet 2005, 37:501506. PubMed Abstract  Publisher Full Text

Wellmer F, AlvesFerreira M, Dubois A, Riechmann J, Meyerowitz E: Genomewide analysis of gene expression during early arabidopsis flower development.
PLoS Genetics 2006, 2(7):10121024. Publisher Full Text

Lohmann JU, Weigel D: Building beauty: the genetic control of floral patterning.
Dev Cell 2002, 2:135142. PubMed Abstract  Publisher Full Text

Liu Z, Mara C: Regulatory mechanisms for floral homeotic gene expression.
Semin Cell Dev Biol 2010, 21:8086. PubMed Abstract  Publisher Full Text

Kaufmann K, Melzer R, Theissen G: MIKCtype MADSdomain proteins: structural modularity, protein interactions and network evolution in land plants.
Gene 2005, 347:183198. PubMed Abstract  Publisher Full Text

de Folter S, Angenent GC: Trans meets cis in MADS science.
Trends Plant Sci 2006, 11:224231. PubMed Abstract  Publisher Full Text

Kaufmann K, Muio JM, Jauregui R, Airoldi CA, Smaczniak C, Krajewski P, Angenent GC: Target genes of the MADS transcription factor SEPALLATA3: integration of developmental and hormonal pathways in the Arabidopsis flower.
PLoS Biol 2009, 7:e1000090. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Riechmann JL, Krizek BA, Meyerowitz EM: Dimerization specificity of Arabidopsis MADS domain homeotic proteins APETALA1, APETALA3, PISTILLATA, and AGAMOUS.
Proc Natl Acad Sci USA 1996, 93:47934798. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Sablowski R: Genes and functions controlled by floral organ identity genes.
Semin Cell Dev Biol 2010, 21:9499. PubMed Abstract  Publisher Full Text

Pellegrini L, Tan S, Richmond TJ: Structure of serum response factor core bound to DNA.
Nature 1995, 376:490498. PubMed Abstract  Publisher Full Text

Honma T, Goto K: Complexes of MADSbox proteins are sufficient to convert leaves into floral organs.
Nature 2001, 409:525529. PubMed Abstract  Publisher Full Text

Theissen G, Saedler H: Floral quartets.
Nature 2001, 49:469471. Publisher Full Text

Immink RG, Tonaco IA, de Folter S, Shchennikova A, van Dijk AD, BusscherLange J, Borst JW, Angenent GC: SEPALLATA3: the 'glue' for MADS box transcription factor complex formation.
Genome Biol 2009, 10:R24. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Immink RG, Kaufmann K, Angenent GC: The 'ABC' of MADS domain protein behaviour and interactions.
Semin Cell Dev Biol 2010, 21:8793. PubMed Abstract  Publisher Full Text

GomezMena C, de Folter S, Costa MM, Angenent GC, Sablowski R: Transcriptional program controlled by the floral homeotic gene AGAMOUS during early organogenesis.
Development 2005, 132:429438. PubMed Abstract  Publisher Full Text

Alon U: An introduction to systems biology. Design principles of biological circuits. Chapman & Hall/CRC; 2006.

Jack T, Fox GL, Meyerowitz EM: Arabidopsis homeotic gene APETALA3 ectopic expression: transcriptional and posttranscriptional regulation determine floral organ identity.
Cell 1994, 76:703716. PubMed Abstract  Publisher Full Text

Zachgo S, Silva EdeA, Motte P, Trbner W, Saedler H, SchwarzSommer Z: Functional analysis of the Antirrhinum floral homeotic DEFICIENS gene in vivo and in vitro by using a temperaturesensitive mutant.
Development 1995, 121:28612875. PubMed Abstract  Publisher Full Text

Kwiatkowska D: Flower primordium formation at the Arabidopsis shoot apex: quantitative analysis of surface geometry and growth.
J Exp Bot 2006, 57:571580. PubMed Abstract  Publisher Full Text

Parcy F, Nilsson O, Busch MA, Lee I, Weigel D: A genetic framework for floral patterning.
Nature 1998, 395:561566. PubMed Abstract  Publisher Full Text

Lenhard M, Bohnert A, Jrgens G, Laux T: Termination of stem cell maintenance in Arabidopsis floral meristems by interactions between WUSCHEL and AGAMOUS.
Cell 2001, 105:805814. PubMed Abstract  Publisher Full Text

Lohmann JU, Hong RL, Hobe M, Busch MA, Parcy F, Simon R, Weigel D: A molecular link between stem cell regulation and floral patterning in Arabidopsis.
Cell 2001, 105:793803. PubMed Abstract  Publisher Full Text

Klipp E, Kowald A, Wierling C, Lehrach H: Systems biology in practice: concepts, implementation and application. WileyVCH; 2005.

Cavelier G, Anastassiou D: DataBased Model and Parameter Evaluation in Dynamic Transcriptional Regulatory Networks.
Proteins 2004, 55:339350. PubMed Abstract  Publisher Full Text

Bundschuh R, Hayot F, Jayaprakash C: Fluctuations and Slow Variables in Genetic Networks.
Biophysical Journal 2003, 84(3):16061615. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Buchler NE, Louis M: Molecular titration and ultrasensitivity in regulatory networks.
Journal of Molecular Biology 2008, 384:11061119. PubMed Abstract  Publisher Full Text

Wigge PA, Kim MC, Jaeger KE, Busch W, Schmid M, Lohmann JU, Weigel D: Integration of spatial and temporal information during floral induction in Arabidopsis.
Science 2005, 309:10561059. PubMed Abstract  Publisher Full Text

Urbanus S, de Folter S, Shchennikova A, Kaufmann K, Immink R, Angenent G: In planta localisation patterns of MADS domain proteins during floral development in Arabidopsis thaliana.
BMC Plant Biology 2009., 9(5) PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gennemark P, Wedelin D: Efficient algorithms for ordinary differential equation model identification of biological systems.

Kimura S, Ide K, Kashihara A, Kano M, Hatakeyama M, Masui R, Nakagawa N, Yokoyama S, Kuramitsu S, Konagaya A: Inference of Ssystem models of genetic networks using a cooperative coevolutionary algorithm.
Bioinformatics 2005, 21:11541163. PubMed Abstract  Publisher Full Text

Kimura S, Hatakeyama M, Konagaya A: Inference of Ssystem models of genetic networks from noisy timeseries data.
ChemBio Inform J 2004, 4:114. Publisher Full Text

Maki Y, Ueda T, Okamoto M, Uematsu N, Inamura Y, Eguchi Y: Inference of genetic network using the expression profile time course data of mouse P19 cells.

Coleman T, Li Y: On the Convergence of Reflective Newton Methods for LargeScale Nonlinear Minimization Subject to Bounds.
Mathematical Programming 1994, 67(2):189224. Publisher Full Text

Coleman TF, Li Y: An Interior, Trust Region Approach for Nonlinear Minimization Subject to Bounds.
SIAM Journal on Optimization 1996, 6:418445. Publisher Full Text

Dormand J, Prince P: A family of embedded RungeKutta formulae.
Journal of Computational Applied Mathematics 1980, 6:1926. Publisher Full Text

Jack T, Brockman LL, Meyerowitz EM: The homeotic gene APETALA3 of Arabidopsis thaliana encodes a MADS box and is expressed in petals and stamens.
Cell 1992, 68:683697. PubMed Abstract  Publisher Full Text

Mizukami Y, Huang H, Tudor M, Hu Y, Ma H: Functional domains of the floral regulator AGAMOUS: characterization of the DNA binding domain and analysis of dominant negative mutations.
Plant Cell 1996, 8:831845. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yanofsky MF, Ma H, Bowman JL, Drews GN, Feldmann KA, Meyerowitz EM: The protein encoded by the Arabidopsis homeotic gene agamous resembles transcription factors.
Nature 1990, 346:3539. PubMed Abstract  Publisher Full Text

van Dijk AD, ter Braak CJ, Immink RG, Angenent GC, van Ham RC: Predicting and understanding transcription factor interactions based on sequence level determinants of combinatorial control.
Bioinformatics 2008, 24:2633. PubMed Abstract  Publisher Full Text

EgeaCortines M, Saedler H, Sommer H: Ternary complex formation between the MADSbox proteins SQUAMOSA, DEFICIENS and GLOBOSA is involved in the control of floral architecture in Antirrhinum majus.
EMBO J 1999, 18:53705379. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Melzer R, Theissen G: Reconstitution of 'floral quartets' in vitro involving class B and class E floral homeotic proteins.
Nucleic Acids Res 2009, 37:27232736. PubMed Abstract  Publisher Full Text  PubMed Central Full Text