Continuous-time modeling of cell fate determination in Arabidopsis flowers.

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, continuous-time 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 MADS-box 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 well-known 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 [2-7]. Despite these modifications, however, the ABCDE-model 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, [8-10]. Also, a stochastic type of boolean network, and a differential equations model, that can be considered as a first approximation of kinetic-reaction 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 B-type 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 toggle-switch gene network [18]. This type of model allows a quantitative, continuous-time 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 MADS-box 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 [25-28].

The consensus MADS domain protein target site (CArG-box) 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 higher-order 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 [30-33]. 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 protein-interaction 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 SEPALLATA1-SEPALLATA4 (SEP1, SEP2, SEP3 and SEP4). SEP1-3 are expressed in whorl 2-4, 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 2-4, and SEP1-SEP4 in whorl 1-4.

Figure 1

The ABCDE model for flower organ determination in . The figure has to be read column-wise. E.g., in the sepal-whorl 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).

Protein	Activator	Repressor
AP1	[AP1 SEP]	[AG AG]
AP3	[AP1 SEP], [AG SEP], [AP3 PI]
PI	[AP1 SEP], [AG SEP], [AP3 PI]
AG	[AG SEP], [AG AG]	[AP1 AP1]
SHP	[AG SEP]	[PI AP3]
SEP	[AG SEP], [AP1 SEP], [SEP SEP]

Figure 2

Graphical representation of the interactions in Table 1.

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.

Whorl day	2	3	4	5
1 (sepals)	AP1	AP1	AP1	AP1
2 (petals)	AP1,AP3,PI	AP1,AP3,PI	AP1,AP3,PI	AP1,AP3,PI
3 (stamens)	AP3,PI,AG	AP3,PI,AG	AP3,PI,AG	AP3,PI,AG
4 (carpels)	AG	AG	AG,SHP	AG,SHP

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 whorl-specific triggering mechanisms, representing the timed initial activation of the MADS domain proteins by non-MADS factors, such as LEAFY (LFY), and WUSCHEL (WUS) [22,23]. (b) After activation, MADS box gene transcription is regulated by auto-regulatory 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 non-functional 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 xthe concentration of monomer i and by [xx] the concentration of the dimer of proteins i and j, we have the following balance equation

The proteins represented by the variable xare listed in Table 3.

Table 3

Notation for the variables that represent the MADS proteins.

variable	MADS
x1	AP1
x2	AP3
x3	PI
x4	AG
x5	SHP
x6	SEP

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 Michaelis-Menten 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, p2(t, w) and p4(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 time-and whorl-dependent 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 Michaelis-Menten 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 quasi-steady 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 γ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 whorl-specific 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 [102, 103] nM, and in [44] values used for Km (there named )) are in the range of [10, 102] 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∈ [10-3, 1] nM-1min-1 and K∈ [10-3, 10] min-1[45]. The relative interaction strengths between dimerising proteins are based on expert knowledge:

where {xx} denotes the value of K/Kcorresponding to the dimer [xx]. Based on this information, we fix the values of these parameters at K= 1, γ1 = 1, γ2 = 10-1, γ3 = 10-1, γ4 = 10-2, γ5 = 10-2, γ6 = 5.10-3, and γ7 = 5.10-3nM-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 SEP1-SEP4 in whorl 2-4, and for SEP4 in whorl 1. The intensities are averages over the whole meristem. Since we need whorl-specific 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 103nM (in [45] a range of 102-104 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·104μm3, the petals of 2.7·104μm3, the stamens of 2.9·104μm3, and the carpels of 1.1·104μm3.

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(t) the organ volume. Further,

with the percentage of expression (1% when there is no expression, 100% when the gene is expressed), and α(t) a scaling factor. To determine α(t), equation (8) is inserted into (7), which yields the expression

Network decoupling

With the γ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 [48-51]. 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 x2, .., x6 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 x2(t), ..., x6(t) and are substituted in (10) so that in the resulting equation x1 is the only variable. This equation is integrated and by fitting the calculated values of x1 to the data for x1, 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, 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 quasi-steady 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(in days), and x,(t, 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 gradient-based 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 Dormand-Prince 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 Kwere chosen such that for typical concentration values the Michaelis-Menten 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.

Results

The estimated parameter values are listed in Table 4.

Table 4

Identified parameters for model (5). βis the maximal transcription rate in nM day-1 for the j th Michaelis-Menten function of gene i, Kmthe corresponding half-maximal activation in nM, dthe decay in day-1, and Pthe trigger in nM day-1.

β1,1	6.6e4	Km1,1	10	Km6,1	5.7e2
β2,1	3.3e4	Km1,2	3.7e2	Km6,2	20
β2,2	1.2e2	Km2,1	6.1e2	Km6,3	47
β2,3	1.2e2	Km2,2	1.1e3	dc1	71
β3,1	1.5e3	Km2,3	1.1e2	dc2	3
β3,2	38	Km3,1	3.1e2	dc3	48
β3,3	38	Km3,2	6.3e2	dc4	5e2
β4,1	8.8e3	Km3,3	63	dc5	4
β4,2	1.5e4	Km4,1	10	dc6	16
β5,1	4.1e2	Km4,2	1e2	P2	4.5e4
β6,1	2.3e3	Km4,3	1.1e3	P4	3.4e3
β6,2	44	Km5,1	1e3
β6,3	6e3	Km5,2	10

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·103, 7·104] 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, 103] nM-1, with the exception of Km2,2 and Km4,3, which are only slightly higher. All decay rates are within the range of [1, 1.5·102] day-1, except dc4, 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 2-5, but still a factor 5 lower than the on-level. 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 log-scale visually magnifies errors in the low concentration regime. Since the assumption that the off-genes have 1% of the concentration of the on-genes 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 [x1, .., x6].

where x is the simulated and xthe measured gene expression. In (13) we take on purpose not the squares of the differences, since epsilon should measure whether xis systematically over-or underestimated over five days. Moreover, the differences x-xhave been normalized by xin 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 on-level. AG has a good fit for the on-levels 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 0-3, and the low levels at days 4-5. 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 B-gene PI become so small that they are not visible in Figure 6, due to the logarithmic scale. In the stamens the B-genes are off, and since in this model SHP is repressed by the AP3-PI 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 wild-type expression levels, and they are shown to compare the mutant dynamics with.

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 1-4 will have sepals-petals-petals-sepals 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 yeast-two-hybrid 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 γ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.

Removed dimer	Phenotypic alteration
[AG AG]	stamens → petals, carpels → sepals
[AG SEP]	stamens → petals, carpels → sepals
[SHP SEP]	carpels: no ovules
[AP1 SEP]	sepals → carpels, petals → stamens
[AP3 PI]	petals → sepals, stamens → carpels
[AP1 AP1]	sepals → carpels, petals → stamens
[SEP SEP]	no development of floral organs

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 high-quality 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 well-defined 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 'wet-lab' experiments and 'dry-lab' 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 co-wrote the article. All authors read and approved the manuscript.

Additional file 1

Sensitivity analysis of the optimal parameters with respect to binding affinities.

Click here for file

Additional file 2

Discussion on the simulation results of the mutants.

Click here for file