### 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.

This study is based on two commonly used rate law formulations: the *Hill* rate law [2, 12], defined by Eq. 1, and the *artificial neural network* (*ANN*) rate law [13], defined by Eq. 2.

\begin{array}{c}\frac{d{x}_{i}}{dt}={\widehat{\alpha}}_{i}{\displaystyle \sum _{j}^{n}}\left|{\omega}_{ij}\right|{r}_{i}\left({x}_{i}\right)-{\beta}_{i}{x}_{i},\phantom{\rule{2.36043pt}{0ex}}with\\ {r}_{i}\left({x}_{j}\right)=\left\{\begin{array}{cc}{{x}_{j}}^{{n}_{ij}}/\left({{x}_{j}}^{{n}_{ij}}+{{\gamma}_{i}}^{{n}_{ij}}\right)\hfill & if{\omega}_{ij}>0\hfill \\ 1/(1+({x}_{j}/{\gamma}_{i}){n}_{ij}\hfill & if{\omega}_{ij}<0\hfill \end{array}\right.\end{array}

(1)

\frac{d{x}_{i}}{dt}={\widehat{\alpha}}_{i}\frac{1}{1+\mathsf{\text{exp}}\left({\gamma}_{i}-\sum _{j}^{n}{\omega}_{ij}{x}_{j}\right)}-{\beta}_{i}{x}_{i}

(2)

**-** *x*_{
i
}*, x*_{
j
} ∈ {*x*_{1}*, x*_{2}*, ..., 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 {\mathbb{R}}^{+}:*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|\phantom{\rule{2.36043pt}{0ex}}\in \phantom{\rule{2.36043pt}{0ex}}{\mathbb{R}}_{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*_{1}*, x*_{2}*, x*_{3}, 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.