ECB-ART-41248BMC Syst Biol 2009 Aug 23;3:83. doi: 10.1186/1752-0509-3-83.
Show Gene links Show Anatomy links
Monte Carlo analysis of an ODE Model of the Sea Urchin Endomesoderm Network.
BACKGROUND: Gene Regulatory Networks (GRNs) control the differentiation, specification and function of cells at the genomic level. The levels of interactions within large GRNs are of enormous depth and complexity. Details about many GRNs are emerging, but in most cases it is unknown to what extent they control a given process, i.e. the grade of completeness is uncertain. This uncertainty stems from limited experimental data, which is the main bottleneck for creating detailed dynamical models of cellular processes. Parameter estimation for each node is often infeasible for very large GRNs. We propose a method, based on random parameter estimations through Monte-Carlo simulations to measure completeness grades of GRNs. RESULTS: We developed a heuristic to assess the completeness of large GRNs, using ODE simulations under different conditions and randomly sampled parameter sets to detect parameter-invariant effects of perturbations. To test this heuristic, we constructed the first ODE model of the whole sea urchin endomesoderm GRN, one of the best studied large GRNs. We find that nearly 48% of the parameter-invariant effects correspond with experimental data, which is 65% of the expected optimal agreement obtained from a submodel for which kinetic parameters were estimated and used for simulations. Randomized versions of the model reproduce only 23.5% of the experimental data. CONCLUSION: The method described in this paper enables an evaluation of network topologies of GRNs without requiring any parameter values. The benefit of this method is exemplified in the first mathematical analysis of the complete Endomesoderm Network Model. The predictions we provide deliver candidate nodes in the network that are likely to be erroneous or miss unknown connections, which may need additional experiments to improve the network topology. This mathematical model can serve as a scaffold for detailed and more realistic models. We propose that our method can be used to assess a completeness grade of any GRN. This could be especially useful for GRNs involved in human diseases, where often the amount of connectivity is unknown and/or many genes/interactions are missing.
PubMed ID: 19698179
PMC ID: PMC2739852
Article link: BMC Syst Biol
Species referenced: Echinodermata
Genes referenced: alx1 arid3a ets1 gata6 hnf6 LOC100887844 LOC100893907 LOC115925415 LOC575170 LOC583082 LOC592057 LOC594353 onecut2 otx2 pmar1 tbr1
Article Images: [+] show captions
|Figure 1. Topology of the Endomesoderm Network. The Endomesoderm GRN, as used in this analysis. Reproduced from  (version of December 2007) with permission from E.H. Davidson. Horizontal lines with bent outward arrow indicate genes. Arrows towards genes indicate activators of transcription, barred lines incident to genes indicate inhibitors of expression. Each gene is assumed to be exclusively regulated by the interactions shown. Protein interactions are indicated by circles. For details of the model, please refer to Additional Files 1 and 6. Genes or proteins without any input shown were set manually in simulations. Cis-regulatory inputs to genes that have not been identified by experimental analysis are denoted 'Ubiq', 'Mat' denotes amternal inputs, 'nucl.' indicates nuclearization and 'X' refers to the source of β-catenin. The different territories simulated contain the same set of possible interactions but due to temporally differing inputs, different subsets are realized in each territory.|
|Figure 2. Schematic summary of the detection of TPEs. The Endomesoderm GRN was derived from experimental data (step 1), from this topology we derive an ODE model (2) which is then repeatedly simulated with randomly sampled parameter sets under different perturbations (3). Exemplary, the possible effect of knock down of an activator (a-KD) and an inhibitor (i-KD) on a genes expression are shown. Qualitative topological perturbation effects, the sign of the difference between control and a-KD resp. i-KD indicated by the red and green arrows, are detected from these simulations (4) and these effects compared to experimental data (5). For the comparison, quantitative experimental data needs to be transformed to qualitative data (6). In order to rank the agreement between network topology and experimental data, the same method is applied to randomized networks (R). Alternatively, the resulting agreement can be compared to the agreement expected given that the correct topology is used, as described in section 'Inference of the Reliability of the Proposed Heuristic' and sketched in Figure 3.|
|Figure 3. Validation of our Method. Schematic representation of the approach we use to evaluate the heuristic we propose. The topology of a given GRN is inferred from in-vivo data (red arrow). To judge this step, we propose our heuristic as described in text and Figure 2 (yellow arrow). The comparison of in-vivo data and detected TPEs yields some measure of matching effects (left light blue arrow). In order to evaluate this, we simulate in-silico data that perfectly matches an underlying topology (green arrow). Comparison between in-silico data and detected TPEs indicates the reliability of the heuristic (lower light blue arrow).|
|Figure 4. Construction of ODEs from network topology. Schematic representation of the conversion of a GRN topology to an ODE model. In the diagram to the left, A, B, C and X denote genes. In the Boolean formula in the center, the variables represent the gene products, translation and transcription are combined to one step. In the differential equation to the right, X indicates concentration of mRNA of gene X, while A, B and C are protein concentrations, small letters are kinetic parameters. The topology is used to derive a Boolean formulation of the regulatory inputs to a gene. This Boolean formulation is automatically translated to an ODE model. For details, see section Methods.|
|Figure 5. Simulation of differential expression. Timecourse simulation for Alx1 (left) and Otx (right) expression in different territories up to 70 hours post fertilization (hpf) in arbitrary units (au). Parameters of the ODE model are not fitted to reproduce experimental data. Spatial expression patterns are nevertheless generated in our model: Expression in the total embryo (blue) is the sum of expression in endoderm (purple), mesoderm (brown) and PMC (green) expression. The plotted values are the means of 800 simulations using randomly sampled parameter sets (error bars not shown). In the left part, total Alx1 mRNA abundance equals abundance in PMC; abundance in endoderm and mesoderm are 0. The simulation results show differential expression in the different territories.|
|Figure 6. Effects of perturbations in Pmar1 expression show effects of double-negative regulation of TBr, Alx1 and Ets1. Simulated timecourse of HesC, TBr (left), and Alx1, Ets1 (right) abundance under Pmar1-MOE, Pmar1-KD and control conditions in arbitrary units up to 70 hours post fertilization. The values plotted are the means of 800 simulation results using different parameter sets (error bars not shown). HesC, the only direct target of Pmar1 is inhibited by Pmar1. The inhibitory role of HesC on TBr, Ets1 and Alx1 causes these genes to react as if activated by Pmar1.|
|Figure 7. Scatterplots illustrate the effect of Pmar1-perturbations. Scatterplots of HesC (left) and Alx1 (right) mRNA abundance. Simulated abundances at timepoint t = 25 hpf for all 800 parameter sets are plotted. Red dots represent simulation results under unperturbed (x-axis) and Pmar1-KD (y-axis) conditions. Green dots correspond to simulation results under unperturbed (x-axis) and Pmar1-MOE (y-axis) conditions. In these examples, the different perturbations are clearly distinguishable. Due to the double-negative regulation (Pmar1 inhibits HesC inhibits Alx1), the perturbations in Pmar1 cause converse effects in the expression of HesC and Alx1.|
|Figure 8. Territory-specific perturbation effects. Topological perturbation effects of different genes for Alx1-KD, Dri-KD, Ets1-KD, Hnf6-KD and Pmar1-KD at different timepoints and for different embryonic territories (all territories, endoderm, mesoderm, PMC). Rows show the effects of the indicated perturbation on the expression of the respective gene at time points 14, 19, 25, 33, 45 and 66 hpf. For each timepoint, the 4 columns indicate the effects for (from left to right) 1.) the combination of endoderm, mesoderm and PMC, 2.) endoderm alone, 3.) mesoderm alone, 4.) PMC alone. Red indicates a decrease in expression, green an increase in expression, white indicates no effect of the perturbation on the respective gene. Gray fields indicate genes that are not expressed in the given territory.|
|Figure 9. TPEs discriminate territory-specific downstream differentiation genes. Groups of genes driven by the Endomesoderm GRN (genes from the lower part of Figure 1) are differently effected by upstream perturbations in simulations. These effects are consistent with the embryonic territory these genes are associated with. The top 6 genes are associated with PMC, the next 3 with mesoderm and the 2 bottom genes with endoderm. TPEs for different timepoints were collated. Red indicates decrease in expression, green indicates increase in expression, dark green indicates a late increase in expression, white indicates no change in expression in response to perturbation.|
|Figure 10. Negative regulation of GataE by Hox increases amount of TPEs of Hox perturbations. Topological perturbation effects of GataE and Hox knockdown on expression of different genes as detected in simulations. Green indicates an increase in expression, red indicates a decrease in expression, white indicates no change in expression. Notice that most genes, including GataE are effected by both GataE-and Hox-KD but that Hox expression is not effected by GataE-KD.|
|Figure 11. Detailed comparison of TPEs and experimental data. Matches and mismatches between experimental data and TPEs. The matches and mismatches for the different timepoints were combined by assigning +1 to matches,-1 to mismatches and calculating the mean for each perturbation/gene pair. Means close to 1 are green, those close to -1 are colored red. Means close to 0 (equal matches/mismatches) are white. Gray cells indicate perturbation/gene pairs for which no experimental data was available. Expression changes for the perturbed gene were not included in the analysis, because translation was blocked in experiments and transcription was blocked in simulations. For updated experimental data, see .|
References [+] :
Ben-Tabou de-Leon, Modeling the dynamics of transcriptional gene regulatory networks for animal development. 2009, Pubmed