Evaluating a common semi-mechanistic mathematical model of gene-regulatory networks

Rate law

The main assumption behind automated GRN model inference from timecourse gene expression data is that such data contains sufficient information to generate models that capture the essential mechanistic characteristics of the underlying biological GRN system. A common strategy for modeling and simulating dynamic GRNs is based on nonlinear ordinary differential equations (ODEs) that are derived from standard mass-balance kinetic rate laws [2]. The ODEs in a GRN model relate changes in gene transcripts concentration to each other (and possibly to an external perturbations). Such models consist of one ODE for each gene in the GRN, where each equation describes the transcription rate of the gene as a function of the other genes (and of the external perturbations). The parameters of the equations have to be inferred from the expression time-course data. ODE GRN models are similar to metabolic models that are formulated based on enzyme kinetics, where each rate law approximates a series of elementary chemical steps. Here, the rate laws are one level of complexity above that and represent a series of enzymatic steps. Because these rate laws combine mechanistic details into a small set of model parameters, they are sometimes referred to as "lumped" or "semimechanistic" models. In a sense, these models are neither fully mechanistic nor purely phenomenological.

- x _i , x _j ∈ {x₁, x₂, ..., x _n }: Time-dependent transcript concentration of gene i and j, respectively, where n is the total number of genes in the GRN system;

- dx _i /dt ∈ ℝ: Total rate of x _i change at time t;

- ${\widehat{\alpha }}_i\in {ℝ}^{+}:$Maximal synthesis rate of transcript x _i ;

- ω _ij ∈ ℝ: Type of synthesis regulation of transcript x _i by x _j , such that

ω _ij >0: Synthesis activation of x _i by x _j ;

ω _ij <0: Synthesis repression of x _i by x _j ;

ω _ij = 0: No synthesis regulation of x _i by x _j .

-$\left|{\omega }_{ij}\right|\in {ℝ}_0^{+}:$Relative weight of synthesis-regulatory influence of x _j on x _i ;

- γ _i : Activation/repression coefficient of gene i; γ _i ∈ ℝ for ANN, and

γ _i ∈ ℝ⁺ for Hill;

- n _ij ∈ ℝ⁺: Hill coefficient controlling the steepness of the sigmoidal regulation function; and

- β _i ∈ ℝ⁺: Degradation rate constant modulating the degradation rate of x _i .

Both rate laws have been shown to represent essential characteristics of biological processes [2, 8, 12–15]. They capture a maximal synthesis rate (${\widehat{\alpha }}_i$), sigmoidal (saturable) kinetics, and an activation/repression threshold (γ _i ). For nij <1, the Hill rate law represents Michaelis-Menten kinetics. The rate law versions shown in Eqs. 1 and 2 assume additive input processing and a linear transcript degradation rate (β _i x _i ) that depends only on the concentration of the target gene's product. These assumptions are not set in stone; the rate laws may be adapted to capture multiplicative input processing and a non-linear degradation rate which may depend on various influences. Variations that capture basal transcript synthesis and input delays are also possible [8, 2].

Like in other comparable GRN rate laws (e.g. the synergistic-system [16]), the omega parameter (ω _ij ) represents two distinct biological concepts simultaneously; a discrete as well as a continuous concept. On one hand, it defines the nature or type of synthesis regulation between two genes i and j: if ω _ij >0, then gene j activates synthesis of transcript x _i , if ω _ij <0, then gene j represses x _i synthesis, and if ω _ij = 0, then gene j does not regulate transcript x _i at all. Hence, the totality of all ω _ij parameters determines the transcript synthesis-regulatory network structure of the GRN. On the other hand, the quantity |ω _ij | defines the strength or influence of a regulator gene j on its target or effector gene i. When we employ automated reverse-engineering of GRN models from time-course gene expression data with algorithms like the one illustrated in Algorithm 1, the dual role of ω _ij and its discrete-continuous interpretation has important consequences.

First, because ω _ij needs to be coded as a real number (ω _ij ∈ ℝ), the chances of a typical parameter estimation algorithm to find ω _ij = 0 are very small. Thus, reverse-engineering algorithms like the one discussed below have a tendency to infer only non-zero values for |ω _ij |, representing fully connected network structures. Fully connected GRN network structures are at best very difficult to interpret biologically, at worst meaningless.

Second, because typical GRN model formalisms like the ones in Eqs. 1 and 2 contain additional parameters to represent other quantitative aspects of GRN systems, reverse-engineering algorithms tend to "balance" the quantitative values of all parameters. This means that only the relative quantities |ω _ij | are important, not their absolute values! It is important to understand this property, as this lies at the heart of this study.

Third, in the absence of a clear understanding of the effect ω _ij has in the inference process, there is a danger that modelers choose large omega intervals in their algorithms. This, of course, adds additional computational load because it increases the size of the search or solution space.

Inference algorithm

Once one has chosen a rate law or model formalism to represent a GRN, one needs to determine the concrete values of the model's parameters - the parameters that describe the network structure, and the parameters that represent other aspects of the modeled GRN system. If these parameters are not known, they are typically inferred by reverse-engineering or parameter estimation algorithms like the one defined by Algorithm 1.

Given stimulus-response gene expression time-course data, D, and a particular model formulation, M, Algorithm 1 determines concrete parameter values. The algorithm iterates over three main steps:

1. An optimizer algorithm that generates candidate model parameter values by attempting to minimize the training error, E.

2. An ODE solver component that numerically integrates the model equations using the initial values of the time series in the training data set, D.

Input: M ← Model equations; L ← Parameter limits; G ← Network topology; D ← Training data; ε ← Error threshold

Output: P ← Parameter values; E ← Training error;

(* Initialize and process: *)

S ← Simulation data (* Initialize *)

E ← ∞ (* Initialize *)

repeat

P ← Optimize(L, E) (* Parameter values *)

S ← SolveODE(M, P, D) (* Solve model *)

E ← Error(S, D) (* Determine error *)

until E < ε ;

Algorithm 1: Basic reverse-engineering algorithm. The network topology, G, is an optional input. In this study, we experiment with known network topology only.

3. A component that computes the simulation error, E, based on the gene expression time-course data in the training data set, D, and the predicted or simulated data, S, determined by the ODE solver.

GRN modeling and simulation software tools implement various features that realize the elements listed above. The diagram in Figure 2 depicts the basic components and "workflow" of a typical GRN model reverse-engineering algorithm. The cloud shape on the left represents the GRN system under study. In this simplified example, the GRN system has 3 genes corresponding to the model variables x₁, x₂, x₃, respectively. The series of dots labeled Experimental Data illustrates the three gene-expression time-series that have been experimentally obtained from the GRN system over 8 consecutive time points. We refer to this as the training dataset. The right part of Figure 2 shows the main reverse-engineering loop. The network shape on the right (labeled GRN Model) is a graphical depiction of a algorithm-generated concrete candidate model with gene-regulatory interaction links between the three genes of the system. The algorithm simulates the system based on the candidate model and the initial condition of the dataset (arrow labeled simulate). The simulation produces a predicted or simulated dataset (curves labeled Predicted Data). The experimental and predicted data are then compared (diamond shape) to assess the quality of the candidate model. If the quality is deemed acceptable, the candidate model is retained as the final candidate model. The final candidate model is still subject to validation on independent data (this is not depicted in the diagram).

The simulation of the predicted time series depicted in Figure 2 involves the numeric integration of the model equations. In terms of computational effort, the ODE solver step accounts for approximately for 80% of the total computing time of Algorithm 1. The reverse-engineering process terminates, when the training error drops below a pre-defined error threshold, or when a maximum number of model evaluations is reached.

For this study we have employed the GRN modeling and simulation tool MultGrain/MAPPER. The tool was developed as a part of the European FP7 project Multiscale Applications on European e-Infrastructures (MAPPER) [17]. The goal of MAPPER was to develop a general framework and technology facilitating the development, deployment and execution of distributed multiscale modeling and simulation applications [18, 19]. Based on tools and services developed in the MAPPER project, MultiGrain/MAPPER realizes the GRN model reverse-engineering process (Figure 2) based on a multi-swarm particle swarm optimization algorithm.

In order to estimate the model parameters, we used the our own implementation of a multi-swarm particle swarm optimization [20] (PSO) algorithm. PSO is a population-based meta-heuristic inspired by the flocking, schooling or swarming behavior of animals. Two main advantages of this method include that it optimizes continuous variables and it has the ability to avoid getting stuck in local minima by using a multi-swarm approach which successively swaps particles across each swarm after a fixed number of iterations in order to increase the "genetic" diversity of the overall swarm. The PSO parameters were set according to the guidelines of Pedersen et al. [21], who performed a meta-analysis of the PSO algorithm, testing its performance for a wide range of parameter values.

Training errors synthetic systems

First, we consider the training errors of the GRN models derived from the synthetic GRN systems in Tables 2 and 3. The training error's means and standard deviations are shown in the bottom-right corner of the tables.

From the average mean training errors, we notice that both sets of Hill models have an average mean training error close to 0.14. This is comparable to average mean training error of the ANN model obtained from the ANN system's data. However, the mean training error (0.3554) of the ANN model obtained from the ANN system's data is more than twice that value. Since the ANN rate law (Eq. 2) has fewer parameters than the Hill rate law (Eq. 1), and hence a smaller degree of freedom, it is harder for the ANN models to fit data obtained from Hill systems. Hill models, on the other hand, can adapt easier to data generated from ANN systems.

While above observations are interesting, the most important information in the context of our investigation relates to the groups of 4 error values for a given system/model combination (e.g. 4 values highlighted in bold font in Table 2), as well as to entire columns of error values (e.g. blocks of errors in italic font in Table 3). Tables 2 and 3 highlight four horizontal groups of 4 training errors in bold; these groups have a standard deviation higher than 0.010. If anything, one would expect the errors to get smaller for larger omega intervals (from left to right), because larger omega intervals relate to a larger solution space. However, in most cases such a pattern is not observed. Indeed, even for the training errors in the bottom three rows in both tables (these were obtained from data of the three systems with large omega values: −20 and +20 for repression and activation, respectively), we cannot find a general improvement of training error for increasing omega intervals. For example, in Table 2 the two horizontal groups of 4 training errors highlighted in italic font do not show a clear pattern of decreasing training errors.

Furthermore, when we look at the profiles of the training errors in the columns of both tables, we notice a good pair-wise similarity of training errors (at least within the four columns relating to the same system/model combination). This is illustrated by two columns highlighted in italic font in Table 3. This means that models inferred with different omega intervals show similar training errors for corresponding data sets. There does not seem to be an advantage of using larger omega intervals in the reverse-engineering process.

Validation errors synthetic systems

The validation error of the inferred models characterizes the predictive power of the models. Tables 4 and 5 show the validation errors of the models inferred from the data of the synthetic systems depicted in Figure 3. The mean validation error of all 192 models inferred from the synthetic systems' data is 0.263 with a standard deviation of 0.127 (not shown in tables). So the mean validation error across all models is ca. 34% higher than the mean training error. The variation of the validation errors is similar to that of the training errors (Tables 2 and 3).

The average mean validation errors are consistent with the averages for the mean training errors, in that, the ANN models' predictive performance on the Hill system's data is much poorer than that of the other three models. In fact, the validation errors reveal that inferring models from data that was obtained from systems that were created with the same rate law (as the model), constitutes a considerable bias. The average mean errors for ANN models obtained from ANN system data, and for Hill models from Hill system data are quite low and similar. However, with mixed configurations (different rate law for system and model), we get much higher average mean validation errors. This relates to the important but frequently ignored issue of the modeling error. The modeling error is due to the fundamental imperfections that arise when we make abstractions of reality in the form of mathematical or computational models. A model, any model, is by definition an approximation of reality [25]. The modeling error quantifies how well the abstraction approximates reality. Conceptualizing a complex phenomena such as GRN systems as a mathematical or computational model is a relatively new modeling abstraction. More research is required to understand how to assess the modeling error in such approaches.

Looking at the data in Tables 4 and 5 in detail, we notice that things are less homogeneous than for training errors. This is to be expected, as predicting the time-courses for unseen stimuli is a much harder task than predicting the timecourses for known inputs. In Table 4 the groups for which the within-group standard deviation is greater than 0.075 are highlighted in bold font. Surprisingly, there are many such groups in Hill/ANN model/system configurations. Still, in terms of the hypothesis we are testing, most groups of four do not show a pattern of decreasing validation error with increasing omega intervals. For example, the two horizontal groups of four validation errors highlighted in italic font in Table 5 illustrate two sets of validation errors that do not vary across the omega interval settings. Indeed, in some cases there is even an increase of error - and in other cases a slight decrease. Likewise, when we look at the vertical validation error profiles in columns (e.g. the two columns highlighted in italic font in Table 5), we notice a general pair-wise similarity for each model group. These observations confirm our hypothesis that the absolute size of the interval for ω _ij is not critical. Even when data is generated with large ω _ij values, the reverse-engineered models can approximate the data equally well with small and large ω _ij ranges.

Training and validation errors yeast system

In terms of the mean training error, the two models perform almost identically. But the mean validation error of the ANN model is nearly twice that of the Hill model! This difference in predictive power is quite remarkable, even though we are testing only four omega conditions. We also observe that the variation (standard deviation) in the ANN model performance (validation error) is much higher than that of the Hill model. Clearly, the Hill rate law has more parameters and hence is more likely to fit complex curves. Still, that the ANN model mean validation error is nearly 100% higher than that of the Hill model (when the mean training errors are similar), seems to be an important observation.

We now analyze how the training and validation performance depends on the omega intervals. We observe essentially a similar pattern as in the evaluation of the synthetic systems. For the two groups of four training errors in Table 6 there seems to be hardly any variation in training error from smaller to larger omega intervals. In the four validation errors of the Hill model, we see a minor variation, but a slight rise in error as we move to larger omega intervals (if anything, the error should become smaller, as more solution possibilities are being explored). And in the validation errors of the ANN model, we notice a considerable variation in validation errors but no pattern of decrease in validation error from smaller to larger omega intervals. So overall, this seems to corroborate the results derived from the synthetic systems in Tables 2 and 3 (training errors), and Tables 4 and 5 (validation errors). It seems, that choosing large (and ad hoc) omega intervals does not make a real difference.