In this section we present examples of the intended use of RMBNToolbox. In the first example we generate a large model for a genetic regulatory network that can be used to produce ground truth data for a microarray simulation [25]. In the second example the structure and stoichiometry of a metabolic network are known, and the kinetic laws are randomly generated. Furthermore, the example demonstrates how metabolic fluxes in a steady state can be decomposed by elementary flux modes. The third example studies network stability using a control theoretic approach. The example generates small networks for which the network structure determines the stability. All the Matlab scripts that are used to generate the following example networks can be found in the examples folder of RMBNToolbox. All the generated example networks can be downloaded as additional files of this article.

### Gene regulatory network

In gene regulatory networks a set of genes produce proteins called transcriptional regulators. Transcriptional regulators bind to the promoter areas of genes, thereby activating or inhibiting their transcription. Most of the genes do not produce transcriptional regulators but their functions may be related to other processes, such as metabolism or cellular growth. Transcriptional regulators are usually thought as the key to the cellular control. In this example we produce a large network model with simple structural characteristics [see Additional file 1]. The model mimics a gene regulatory network.

In the generated network there are 1000 transcription reactions which produce one product each. The total of 200 of the products act as transcriptional regulators which control the network by activating and inhibiting the transcription reactions. Each of the transcription reactions has one activatory and one inhibitory regulator which are selected randomly from the 200 regulators.

The synthesis of the gene products is modeled similarly in all cases. The modeling concentrates on the kinetics of transcription and uses the rate law suggested in [23]. The amount of the protein product, for which the gene is a precursor, is assumed to be equal to the produced transcript. Because the number of gene copies is restricted and only a limited number of regulators are able to bind simultaneously, the kinetic law saturates both with the amount of activator and inhibitor. The rate of transcription is

r={V}_{basal}\frac{{K}_{I}^{{n}_{I}}}{{I}^{{n}_{I}}+{K}_{I}^{{n}_{I}}}\frac{{A}^{{n}_{A}}}{{A}^{{n}_{A}}+{K}_{A}^{{n}_{A}}},

(3)

where *V*_{
basal
}is the rate of transcription in the absence of activators and inhibitors, *I* and *A* represent the concentrations of inhibitor and activator, *K*_{
I
}and *K*_{
A
}represent the concentrations with which the inhibitor and the activator have the effect of half of their maximal effects, and *n*_{
I
}and *n*_{
A
}act as Hill coeffcients. The parameter values are random numbers from the following uniform distributions: *V*_{
basal
}∈ *U* (5,10), *K*_{
I
}∈ *U* (2,3), *K*_{
A
}∈ *U* (1,2), *n*_{
I
}∈ *U* (1,2), *n*_{
A
}

*U* (1,2). The initial concentrations *I* and *A* are random numbers from the uniform distributions *I* ∈ *U* (0, 1) and *A* ∈ *U* (0, 1).

The degradation kinetics of each gene product follows the mass action law

*r* = *k P*, (4)

where *k* is a rate parameter and *P* is the concentration of the gene product. Similarly to the kinetic laws of transcription reactions, the parameter values are unique for each degradation reaction. In this case, the value of parameter *k* is drawn from the uniform distribution *U* (0.01, 0.02).

For a comparison, we additionally simulate a duplicate network which mimics a gene deletion [see Additional file 2]. In the duplicated network, the production of a randomly chosen regulator is stopped by setting the parameter *V*_{
basal
}of its transcription reaction to zero. Otherwise the duplicated network is identical to the original network.

The behavioral differences are illustrated by a time series simulation. After constructed, both networks are exported in SBML format and SBML ODE Solver [42] is used to simulate them. Figure 2 gives an overview to differences in simulation results. For each of the species, the figure shows concentration differences between the two simulations. For each time point *t*, the difference *d* (*t*) is calculated as

*d* (*t*) = *c* (*t*) - *c** (*t*), (5)

where *c* (*t*) and *c** (*t*) are concentration values of the species in the first and in the second simulation, respectively. Although most of the species act similarly in both simulations, there are large and unforeseeable dynamic variations too. The effects of the inactivation of the regulator do not fade away or relax to a constant but the inactivation seems to have complex behavioral consequences.

### Simulation and stoichiometric analysis

In this example, a time series simulation is used to illustrate a result from stoichiometric network analysis. As presented in [43], any feasible steady state flux distribution of a metabolic network is a linear combination of so-called elementary flux modes (EFM). We show using random reaction kinetics that this holds in a small examplary metabolic network model. Extreme pathways [44] and elementary fluxes [45] are similar concepts to EFMs, and they would be equally valid for this analysis.

In metabolic systems, the time derivates of metabolite concentrations **c** can be written as presented in Equation 1, i.e., \frac{dc}{dt}=Sv where *S* is the *m* × *n* stoichiometric matrix with *m* metabolites and *n* reactions, and **v** contains the reaction velocites with *v*_{
i
}≥ 0 for each irreversible reaction *i*. Metabolites are classified to external for which it is assumed that the environment always balances their concentrations **c**_{
ext
}, and internal for which the concentrations **c**_{
int
}are determined by the network. A network is then said to be in a steady state if

\frac{d{c}_{int}}{dt}=0,

(6)

i.e., there is no accumulation or depletion of internal metabolites. Specific reaction velocities (flux distributions) are needed to maintain steady states.

Elementary flux modes describe such reversible and irreversible pathways in the network which maintain steady states when working. In an elementary flux mode, each reaction is assigned with its relative velocity compared to other reactions in the same EFM. EFMs are minimal in the sense that the active reactions in an EFM cannot be a subset of the active reactions in another EFM. Elementary flux modes can be calculated based on a stoichiometric matrix and the respective reaction irreversibilities [43]. Let vector **e** denote an elementary flux mode in which element *e*_{
i
}= 0 if reaction *i* is inactive and *e*_{
i
}≠ 0 if the reaction *i* is active. Further, let the set of all *N* elementary flux modes of the network be in matrix *E* = [**e**_{1}, **e**_{2}, ..., **e**_{
N
}]. Then any flux distribution **v**, which results a steady state into the network, can be described as a linear combination of the EFMs as

**v** = *E β*, *β*_{
j
}≥ 0 if EFM *j* is irreversible (7)

where the vector *β* weigths each of the elementary fluxes by a scalar. The weigths are non-negative for EFMs describing irreversible pathways.

Because the calculation of EFMs uses the steady state assumption, Equation 7 has solutions for *β* only if **v** maintains a steady state. Usually the number of EFMs (columns in *E*) is much larger than the number of reactions (rows in *E*), and therefore unique solutions are rare for Equation 7. However, we can test the existence of solutions by setting up a linear programming problem

\begin{array}{l}\mathrm{max}{1}^{T}\beta \\ \begin{array}{ll}suchthat\hfill & E\beta =v\hfill \\ {\beta}_{j}\ge 0,EFMjisirreversible\hfill \end{array}\end{array}

(8)

The objective function is set to find the maximum of the sum of the weigths. Rather than the maximum value, we are now interested in the existence of any solution. In the following, we utilize the fact that the maximum can be found only if any solutions exist.

A hypothetical example network, illustrated in Figure 3, consists of three species and six reactions [46]. The network structure is provided to RMBNToolbox as a stoichiometric matrix. Initial amounts for metabolites, kinetic rate laws for reactions, and their parameter values are chosen randomly, because we aim at illustrating that an arbitrary steady state flux distribution is a linear combination of elementary flux modes.

Program Metatool [47] is used to calculate the elementary flux modes of the generated network. The network is exported in SBML format [see Additional file 3] and simulated using SBML ODE Solver [42] until a steady state is reached. During the simulation, the flux distribution **v** and time derivates of species amounts \frac{dc}{dt} are sampled for every second. Linear programming problems, as described in Eq. 8, are solved for each flux distribution sample. Figure 4 shows the time derivates and the existence of weights *β* for each sample. After a steady state is reached (i.e., the time derivates of internal metabolites become zero), then the flux distribution is a linear combination of the elementary flux modes (i.e., the weigths *β* are found).

### Network stability

Neither RMBNToolbox nor Systems Biology Markup Language take care of the rationality of the generated network models. Possible unstability of the generated model is a typical issue the user has to consider. In this example we look how to exploit control theory for generating models which are unstable and biologically unreasonable and, on the other hand, stable and biologically more reasonable.

In control theory, a system is stable if it has a bounded response to a bounded input [48]. We concentrate on network models which can be formulated as systems of first-order linear differential equations. Their compact representation form is

in which the vector **c** contains states of variables, and the system matrix *A* determines the system properties. The system is known to be stable if the eigenvalues of *A* have nonpositive real parts, and every eigenvalue with zero as the real part has an associated Jordan block of order one [49, 50].

Next we derive a model for which the user can determine if this stability requirement is fulfilled or not. For this purpose, the biochemical network model described by Equations 1 and 2 needs to be represented in the format of Equation 9.

We begin the model reformulation from kinetic laws, i.e., Equation 2. Because the intended model in Equation 9 is linear, we can utilize for its construction such kinetic laws which are linear too. Kinetic laws of the form first-order mass-action fulfill this need. For example, the kinetic law for reaction *j* is *v*_{
j
}= *k*_{
j
}*c*_{
j
}where *k*_{
j
}is a reaction-specific rate constant and *c*_{
j
}is the concentration of the subsrate. Kinetic laws of this form make it possible to represent the reaction rate vector **v** of Equation 1 by a matrix-vector multiplication

**v** = Γ**c**, (10)

where Γ is a diagonal *n* × *m* matrix storing rate constants *k*_{
j
}, for each reaction *j* = 1, ..., *n*, on its main diagonal.Substituting this to Equation 1, it becomes

\frac{dc}{dt}=Sv=S\Gamma c.

(11)

Multiplication of the time invariant matrices *S* and Γ results to the *m* × *m* matrix *A*. The substition of *A* to Equation 11 brings the network model to the form presented in Equation 9.

During the model construction, matrices *S* and Γ are randomly generated after which they are multiplied to produce the matrix *A*. The eigenvalues of *A* are calculated, and the values of their real parts are examined. The generation is repeated until the eigenvalues indicate the required stable or unstable behavior of the model.

As the first example, we generate an unstable network model [see Additional file 4] by requiring at least one positive eigenvalue for the system matrix *A*. The structure of the generated network is depicted in Figure 5. We note that the network includes two kinds of features that are not reasonable in real biochemical networks and which obviously cause instability for the model. The first unrealistic feature is that the mass balance does not hold: Species S1 is decomposed in two parts, S2 and S3, in reaction R1. Further, reaction R2 converts one S2 molecule back to two S1 molecules. This results that S1 can be decomposed to S2 and S3 without loss of mass and, after each decomposition, the amount of S1 becomes doubled. The second problematic feature is the generation of dead-ends in the model: Species S3 is produced by reaction R1, but it is not consumed by any reaction. Therefore, the amount of S3 increases as long as there is a supply of S1. All this causes unstable behavior for the model, as demonstrated by Figure 6 in which species amounts increase rapidly towards infinity.

In the second example the generated network model [Additional file 5] is stable. We note from Figure 7 that the network does not have the two structural features which appeared in the previous model. Instead, the structure is more treelike and without feedback loops, thereby enabling material flows through the network. Correspondingly, the species amounts show relaxation towards constant values in Figure 8.

## Comments

View archived comments (1)