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.


Result
Based on the analysis proposed in the article, we computed the robustness of expression of different genes against random parameter variations using Monte Carlo simulations.
To this end, we performed multiple simulation runs with different parameter sets, each generating a distinct result. The resulting mRNA concentrations have been averaged for each time point, also recording the variation. Using the coefficient of variation c v = σ/µ, the variation in expression of different genes can be compared. It indicates whether the expression of a gene is highly affected by or robust against parameter changes.
The robustness computed here is not synonymous with robustness as commonly used in biological contexts [1]. In most publications, robustness in a GRN implies that the main features of the system are preserved when small deviations from some natural condition occur (e.g. small changes in individual parameter values). In contrast, we simulate the model using sets of arbitrary parameter values, which might be well out of range of biologically meaningful values.
We computed c v for all genes in the model using 800 parameter sets under unperturbed conditions. For each gene, the highest c v (pertaining to the lowest robustness) of all time points is used as an indicator of the genes robustness. Generally, the robustness is higher at earlier measurement points. This is due to 1 variations in expression of genes upstream of the analyzed genes that take time to propagate to the affected gene. Table 1 gives an overview over the c v of each gene in the network. The c v range from 0.21 to 25.69 for the 800 sampled parameter sets used, indicating that the genes in the network, as it is, differ substantially in their robustness against random parameter changes.
We could neither detect a correlation between node in-degree (the number of regulatory inputs to a particular gene) and robustness nor could we detect a correlation between node out-degree (the number of genes a particular gene regulates) and robustness. The first correlation would indicate that the robustness of expression increases with increased number of regulators, the second would indicate that important regulatory genes are more robust than other genes. We could also not find a correlation between a gene's position in the network (early endomeso, early PMC, late endoderm and mesoderm, late PMC) and its robustness.
The genes of the network that inhibit their own expression, namely SoxB1, GataC, FoxB, TBr, FoxA and Blimp1 are the genes with the lowest c v . This indicates that autoinhibition strongly stabilizes a genes expression. Interestingly, the next five genes in the list (Eve, Apobec, OrCt, GataE, HesC ) all have at least one negative regulatory input.
Besides the robustness of negatively regulated genes, the results indicate that the current architecture does not favor any specific set of genes in their robustness against random parameter changes, that the network is too small to detect significant sets of genes that are expressed with similar robustness or that the parameter values used in the analysis are not in the range in which the system can produce a robust response.

Method
In order to infer the robustness of different components of the Endomesoderm Network model we employed simulation results of the unperturbed model M 0 . Using a number of parameter sets J, we obtained different concentrations s(k, 0, j, t) for each parameter set j in J. For each time point, we computed the mean µ and standard deviation σ 2 assuming that the simulation results are normally distributed around a true concentration value. The relative variation c v = σ/µ provides a comparable measure of variation in relation to the specific mean. We used this relative variation here because expression strength among the genes of the models differs substantially. The relative variation is an indicator of how much parameter values affect the expression of different genes, in other words, how robust the genes are to changes in parameters. To infer an overall robustness of each gene for the entire simulation time, we considered the highest relative variation (lowest robustness) as a lower bound. Table 1. Robustness to random parameter changes.
Coefficient of variation (c v ) as a means of robustness against parameter changes, node indegree (ID) and outdegree (OD) for all genes in the Endomesoderm GRN. Gene