arXiv:1512.05055v1 [] 16 Dec 2015

Inferring Gene Regulatory Network Using
An Evolutionary Multi-Objective Method
Yu Chen

Xiufen Zou

School of Science,
Wuhan University of Technology,
Wuhan, 430070, China
Email: ychen@whut.edu.cn

School of Mathematics
and Statistics, Wuhan University,
Wuhan, 430072, China
Email: xfzou@whu.edu.cn

Abstract—Inference of gene regulatory networks
(GRNs) based on experimental data is a challenging
task in bioinformatics. In this paper, we present a biobjective minimization model (BoMM) for inference
of GRNs, where one objective is the fitting error of
derivatives, and the other is the number of connections in the network. To solve the BoMM efficiently,
we propose a multi-objective evolutionary algorithm
(MOEA), and utilize the separable parameter estimation method (SPEM) decoupling the ordinary
differential equation (ODE) system. Then, the Akaike
Information Criterion (AIC) is employed to select
one inference result from the obtained Pareto set.
Taking the S-system as the investigated GRN model,
our method can properly identify the topologies and
parameter values of benchmark systems. There is no
need to preset problem-dependent parameter values
to obtain appropriate results, and thus, our method
could be applicable to inference of various GRNs
models.

I. I NTRODUCTION
Inference of gene regulatory networks (GRNs) is
important to understand every detail and principle
of biological system, and one of the most popular
methods is to reconstruct GRNs by time-course
data. In the past decades, the Boolean networks
[1], the Bayesian networks [2], [3], and the ordinary differential equation (ODE) networks [4]–
[7], etc., have been proposed to reconstruct GRNs,
and meanwhile, the corresponding algorithms were
also proposed to infer the network topologies and
parameter values.
In this paper, we are devoted to infer the S-system
model of GRN is, a popular nonlinear ODE model

that represents the dynamic process of biochemical
system. It is a nonlinear ODE system
N
N
Y
Y
dXi
h
g
Xj ij ,
Xj ij − βi
= αi
dt
j=1
j=1

i = 1, . . . , N,

(1)
where Xi represents the expression level of gene
i, and N is the total number of genes in the
investigated network. There are totally N ∗(2N +2)
parameters to be specified in an S-system, where
αi , βi ∈ R+ are positive rate constants, and gi,j
and hi,j ∈ R are kinetic constants.
Generally speaking, it is an inverse problem to
infer the S-system using the time course data on
expression levels of genes. A general way is to minimize the normalized errors between experimental
time course data and numerical results [8] , that is,
2
T X
N 
X
Xi,cal (t) − Xi,exp (t)
min err(Θ) =
,
Xi,exp (t)
t=0 i=1
(2)
where Xi,cal (t) is the numerical results of Xi at
time t, and Xi,exp (t) refers to the experimental
value of gene expression level. However, the inverse
problem is usually ill-posed, and an L1 regularizer
could be introduced to penalize the data error [9].
In this way, the minimization model is improved as
min

err(Θ) + λL(Θ).

(3)

Here,
L(Θ) =

N X
N
X
i=1 j=1

|gi,j + hi,j |

(4)

is the L1 norm of parameter vector Θ, and λ is
the penalization parameter. However, it is a challenging task to preset an appropriate value of λ for
GRNs with unknown properties because it is often
problem-dependent.
One of the ways to eliminate the difficulty of
set appropriate values of λ is to convert the singleobjective optimization models to bi-objective models. Liu and Wang [12] proposed a three-objective
optimization model simultaneously minimizing the
concentration error, slope error and interaction error
to find a suitable S-system model structure and
its corresponding parameter values, however, they
transformed it to a single-objective optimization
problem solved by a hybrid differential evolution algorithm, that is, a prior preference was introduced to
the inference process. Spieth et al. [11] took the data
error and connectivity number as two minimization
objective, but the S-system was integrally inferred
because no decoupling method was employed. As
a consequence, the obtained Pareto front contained
numerous Pareto vectors, and it was challenging
to choose the appropriate Pareto solution(s) as the
potential S-system setting(s). Koduru et al. [13] and
Cai et al. [14] simultaneously minimized data error
for several different data sets, but, did not optimize
a metric to make the obtained network sparse.
When GRNs are reconstructed via evolutionary
computation methods, another problem arises with
the heavy evaluation process of candidate parameter
settings. To evaluate a set of parameters generated
by recombination strategies, the ODE system must
be solved via some numerical method, for example,
the Runge-Kutta method. However, all equations
of a coupled ODE system must be simultaneously
solved, and thus, the evaluation process is computationally heavy. To reduce the complexity of
individual evaluation, Tsai and Wang [4] used a
linear Lagrange polynomial to decouple the ODE
system, however, another parameter that has to be
regulated was introduced. Liu et al. [10] developed
a separable parameter estimation method (SPEM)
to decouple the S-System. By these means, each
ODE can be respectively identified by optimization
methods.
In order to eliminate the pruning process for

setting appropriate parameter values, in this paper
we present a multi-objective optimization model for
inference of GRNs, and decoupling the S-system
by the SPEM. Then, a multi-objective evolutionary algorithm is proposed to obtain the network
topology and parameter values. Although SPEM
includes computation process of inverse matrices,
our algorithms can perform well for benchmark
problems, because only small population is needed
to obtain satisfactory results.
The rest of this paper is organized as follows.
Section II introduces SPEM and bi-objective optimization model for GRN inference, and a proposed
multi-objective evolutionary algorithm is introduced
in Section III. Then, Section IV performs a preliminary comparison between the obtained Pareto
sets with results reported in references. Ultimate
inference results obtained by Akaike Information
Criterion (AIC) are characterized in Section V,
and finally, Section VI draws the conclusions and
presents the future work.
II. M ETHOD
A. Decoupling S-systems via the Separable Parameter Estimation Method [10]
Trying to fit the curves of derivatives instead
of the function curves, the SPEM minimize the
difference between dXi /dt and
α̂i

N
Y

gˆ

Xj ij − βˆi

N
Y

hˆ

Xj ij ,

j=1

j=1

where
Θ̂ = (α̂i , β̂i , gˆi,j , hˆi,j )
is a candidate parameter vector of the S-system, and
dXi /dt is approximated by the five-point numerical derivative method. Then, SPEM constructs the
following minimization problem
min Ji (αi , βi , gi , hi ) =



N
N
n
Y
Y
P
h
g
,
Si (tk ) − αi
Xj ij (5)
Xj ij − βi
k=1

j=1

j=1

where gi = (gi,1 , . . . , gi,N ), hi = (hi,1 , . . . , hi,N ),
Si (tk ) is the approximate value of d(Xi )/dt

at time tk . Denoting Γi = (αi , βi ),
[Si (t1 ), . . . , Si (tn )] and
 QN
Q
g
hij
X ij − N
j=1 Xj,t1
 j=1. j,t1
..
..
X̃i = 
.

QN
QN
gij
h
− j=1 Xj,tijn
j=1 Xj,tn

Si



,


=

algorithm (MOEA), each equation is represented by
a combination of bit-string and real variables
x = (bx, rx),

where bx = (bg1 , . . . , bgN , bh1 , . . . , bgN ) and
rx = (g1 , . . . , gN , h1 , . . . , hN ) are respectively bitstrings and real vector of length 2N . If bgi = 1,
we can rewrite (5) as
it means that in the model the parameter gi is not
zero; If bhi = 1, it means that in the model the
min Ji (Γi , gi , hi ) = kSi − X̃i Γi k22 ,
(6) parameter h is not zero. Representing a model via
i
x
can
definitely
express the model topology by bitwhere k · k is the l2 -norm. Note that X̃i does not
strings,
and
consequently,
no threshold is needed to
depend on Γi , Γi can be estimated as
round
a
small
parameter
value
to zero. The MOEA
−1

T
T
to
solve
the
multi-objective
optimization
model (9)
(7)
Γ̂i = X̃i X̃i
X̃i Si .
is described as follows.
Substitute (7) into (6), we can infer each equation
Step1: Randomly generate two populations
of S-system by
[P op, RP op] and [OldP op, OldRP op]
of P opSize individuals. Here P op and
T
T
min Ji (gi , hi ) = SiT [I − X̃i (X̃i X̃i )−1 X̃i ]Si .
OldP op are the binary populations,
and RP op and OldRP op are the real
i = 1, . . . , N (8)
populations.
Evaluate
[P op, RP op]
After the kinetic constants are obtained, correspondand [OldP op, OldRP op] via (9).
ing rate constants can be figure out by (7). Then,
Set [ArchiveP op, ArchiveRP op]
=
the coupled S-system is decoupled, and parameters
[OldP op, OldRP op]
and
can be estimated equation by equation.
POOL
=
RP op;
initialize
gbest = (bgbest, rgbest) to be
B. The Multi-objection Optimization Model
a randomly selected individual in
To infer parameters values equation by equation,
[P op, RP op];
we propose a multi-objective optimization (MO)
Step2: Generate
P opsize
offsprings
to
model
construct the intermediate population
(
Ji (gi , hi )
[IP op, IRP op] and evaluate it.
min
(9)
Step3: Set [OldP op, OldRP op] = [P op, RP op].
kgi k + khi k.
Perform non-dominated sorting on
The first objective minimizes the error for estithe combination of [P op, RP op] and
mation of derivatives, and the second objective is
[IP op, IRP op]. Then, sort it via the
employed to obtain a sparse network, where k · k
dominance rank and values of the second
refers to the l0 -norm of parameter vectors, i.e.,
objective in ascending order;
number of connections in the network. Then, we
Step4: Select P opSize best individuals to uppropose a multi-objective evolutionary algorithm to
date the population [P op, RP op]; upobtain the parameter values.
date [ArchiveP op, ArchiveRP op] and
POOL via [P op, RP op]; randomly select
III. M ULTI - OBJECTIVE E VOLUTIONARY
a non-dominated individual as gbest;
A LGORITHM
Step 5:If the stopping criterion is not satisfied,
Because we employ the SPEM decoupling the
goto Step2; otherwise, output the nonS-system, we can infer ODE equations one by
dominated solutions and the iteration proone. In the proposed multi-objective evolutionary
cess ceases.

Because the SPEM is employed decoupling the Ssystem, there are only 2N parameters to identify
each ODE equation. Note that the second minimization objective as the total number of connections.
There are at most 2N +1 non-domimnated solutions
in the Pareto set. Then, we set the population size
greater than 2N +1, and no diversity mechanisms is
needed to obtain a diverse Pareto front. However,
if we take the norm of parameter vectors as the
second objective that can be any non-negative real
value, to obtain a diverse Pareto front a diversekeeping strategy is necessary. As a consequence, the
time complexity increases. Meanwhile, the obtained
Pareto set could contain several Pareto solutions
with same network topology and similar parameter
values, which does not deserve the increased time
complexity to some extent.
Due to the several discrete values of the second
objective, the population could be absorbed by some
network topology (confirmed by the bit-string bx)
that is easy to locate. Thus, we employ a pool of real
vectors that keep the diversity of obtained solutions
at the early stage of evolution and focus on local
exploitation at the late stage.
A. Generation of Offsprings
To generate candidate solutions, three parents
are randomly selected. According to the method
proposed in [15], three parents are selected as
follows.
• With a probability p1 , two different parents
x1 = (bx1 , rx1 ) and x2 = (bx2 , rx2 )
are selected from the present population
[P op, RP op], and x3 = (bx3 , rx3 ) is selected
from the last population [OldP op, OldRP op];
• otherwise, three different parents x1
=
(bx1 , rx1 ), x2 = (bx2 , rx2 ) and x3 =
(bx3 , rx3 ) are randomly selected from the
Archive [ArchiveP op, ArchiveRP op];
Then, compare x1 with x2 , and initialize the offspring off = (boff , roff ) as the winner.
1) Generation of Binary Offsprings:
To generate the binary offspring boff
=
(bof f (1), . . . , bof f (N )),
we
employ
gbest = (bgbest, rgbest) guiding the search.
∀i ∈ {1, 2, . . . , 2N },

If bx2 (i) = bx3 (i) and bx1 (i) 6= gbest(i) hold,
bof f (i) is set to be gbest(i) with probability
1 − p2 ;
• If bx2 (i) = bx3 (i) and bx1 (i) = gbest(i) hold,
bof f (i) is randomly generated with probability p2 ;
• otherwise, bof f (i) is randomly generated with
probability p2 .
2) Generation of Real Offsprings: The real offspring roff = (rof f (1), . . . , rof f (N )) is generated as follows.
• With probability p3 ,
•

rof f (i) = rx1 (i) + F · (rx2 (i) − rx3 (i));
•

otherwise,
rof f (i) = rof f (i) + F · (pool2 (i) − pool3 (i));

Here pool1 = (pool1 (1), . . . , pool1 (N )) and
pool2 = (pool2 (1), . . . , pool2 (N )) are two real
vector randomly selected from the real vector set
POOL;
B. Update of the Archive
At the beginning, the archive is randomly initialized. During the evolving process, arcw , the worst
archive member whose second objective value is
greatest will be replaced once an offspring off is
generated dominating arcw .
C. Update of the Real Pool
To enhance the global exploration and local exploitation abilities, we employ a real pool POOL
with its corresponding objective vector set Fpool
in the proposed MOEA. POOL is initialized as
the real population RP op at the beginning, and
during the evolving process, a randomly selected
real vector pool ∈ POOL is updated by roff if
it dominates pool via the corresponding objective
vectors.
IV. P RELIMINARY R ESULTS FOR I NFERENCE OF
B ENCHMARK S-S YSTEMS
A. Inference Results of Benchmark S-Systems Obtained via Sole Data Set
To evaluate the performance of our proposed
method, we utilize it inferring four benchmark Ssystems. For each benchmark S-system, we employ

TABLE I: Parameter Settings in the Proposed
MOEA.
Parameter setting
p1 = 4 × (I − 0.5)2
p2 = 0.2
p3 = exp (−I)
F = (I + 0.1) × rand(−1, 1)

3

2.5
Expression Level

P arameter
p1
p2
p3
F

3.5

2
X1

1.5

X

2

X

1

3

0.5

0

where Gen is the generation counter, M axGen is
the maximum generation for the stopping criterion.
We perform 20 independent runs for each ODE
equation of four benchmark systems, and the best
obtained results with smallest data errors are respectively included in Tabs. II, III, IV and V. Both
precise network parameters and obtained results are
listed in the tables, where the inference results are
enclosed by parentheses. For convenience of comparison, the first number in the parentheses is results
obtained by our method, and the second is what
was reported in references. If the inference result
for a run cannot correctly identify all connections,
we take it as a failure run of the proposed MOEA.
Then, the success rates are also included in the last
columns of the tables.
1) Inferring Results of S1: The first benchmark
system S1 is a 3-D system where the parameters are
listed in Tab.II. For initial conditions X1 (0) = 1.9,
X2 (0) = 0.9 and X3 (0) = 3.0, data are sampled
from time 0 to 5 with stepsize 0.1, and the timecourse curves are illustrated in Fig.1 [10]. Correct
network connections for X2 and X3 are always included in the obtained non-dominated sets, however,
for the first equation, 3 of 20 runs cannot correctly
identify the network topology. Curve of X1 in Fig. 1
demonstrates that derivative of X1 quickly decrease
to zero after around 0.5 time unit, which implies
that only five sample point of derivative are not

0

1

2

3

4

5

Time

Fig. 1: Trajectories of variables in S1.
10
X

1

9

X

8

X3

2

X4
7
Expression Level

a population of size 20 to solve the proposed
multi-objective optimization model (9). After 4000
iterations, the obtained non-dominated solutions and
their objective values are saved. The parameter
settings listed in Tab. I are controlled via a process
indicator


Gen
,
I = 0.9 −
M axGen/10

6
5
4
3
2
1
0

0

1

2

3

4

5

Time

Fig. 2: Trajectories of variables in S2.

zero. Because our method is based on fitting of
derivatives, it cannot perform well for this case.
2) Inferring Results of S2: The second benchmark system is a 4-D system where the parameters
are listed in Tab.III. For initial conditions X1 (0) =
10, X2 (0) = 1,X3 (0) = 2 and X4 (0) = 3,
data are sampled from time 0 to 5 with stepsize
0.1, and the time-course curves are illustrated in
Fig.2 [10]. Compared with results presented in [10],
the precisions of parameter values are a bit lower.
However, we can obtained the correct network
topology with 100% success rate except for the
third equation. Considering that in [10] the initial
parameter setting is generated by normal mutation
performed on the known parameter values, we can
conclude that our method is competitive to the
method presented in [10] because we randomly
generate initial parameters at the beginning.
3) Inferring Results of S3: The third benchmark
system is a 5-D system where the parameters are
listed in Tab.IV. For initial conditions X1 (0) = 10,
X2 (0) = 1,X3 (0) = 2, X4 (0) = 3 and X4 (0) = 4,

TABLE II: Precise and Obtained Values for Parameters of Benchmark System S1.
i
1
2
3

αi
12
(12.00,10.97)
10
(9.80,9.01)
3
(2.56,2.94)

gi1
0
(0.00,0.00)
0.5
(0.51,0.59)
0
(0.00,0.00)

gi2
0
(0.00,0.00)
0
(0.00,0.00)
0.75
(0.82,0.76)

gi3
-0.8
(-0.80,-0.97)
0
(0.00,0.00)
0
(0.00,0.00)

βi
10
(9.98,8.81)
3
(2.83,2.17)
5
(4.5,4.92)

hi1
0.5
(0.50,0.60)
0(0)
(0.00,0.00)
0
(0.00,0.00)

hi2
0
(0.00,0.00)
0.75
(0.78,0.89)
0
(0.00,0.00)

hi3
0
(0.00,0.00)
0
(0.00,0.00)
0.5
(0.55,0.50)

success rate
85%
100%
100%

TABLE III: Precise and Obtained Values for Parameters of Benchmark System S2.
i
1
1
3
4

αi
12
(11.80,12.00)
8
(8.00,8.00)
3
(3.00,3.00)
2
(1.90,2.00)

gi1
0
(0.00,0.00)
0.5
(0.50,0.50)
0
(0.00,0.00)
0.5
(0.49,0.50)

gi2
0
(0.00,0.00)
0
(0.00,0.00)
0.75
(0.76,0.75)
0
(0.00,0.00)

gi3
-0.8
(-0.79,-0.80)
0
(0.00,0.00)
0
(0.00,0.00)
0
(0.00,0.00)

gi4
0
(0.00,0.00)
0
(0.00,0.00)
0
(0.00,0.00)
0
(0.00,0.00)

βi
10
(10.00,9.97)
3
(3.00,3.01)
5
(5.10,5.01)
6
(5.80,6.00)

hi1
0.5
(0.50,0.50)
0
(0.00,0.00)
0
(0.00,0.00)
0
(0.00,0.00)

12

hi2
0
(0.00,0.00)
0.75
(0.75,0.75)
0
(0.00,0.00)
0
(0.00,0.00)

hi3
0
(0.00,0.00)
0
(0.00,0.00)
0.5
(0.51,0.50)
0
(0.00,0.00)

hi4
0
(0.00,0.00)
0
(0.00,0.00)
0.2
(0.20,0.20)
0.8
(0.80,0.80)

success rate
100%
100%
90%
100%

1.4
X1
X2

10

1.2

X

3

1

X5

8

Expression Level

Expression Level

X4

6

0.8

0.6
X

1

4

X2

0.4

X3
2

X

4

0.2

X

5

0

0

2

4

6

8

10

Time

0

0

0.1

0.2

0.3

0.4

0.5

Time

Fig. 3: Trajectories of variables in S3.

Fig. 4: Trajectories of variables in S4.

data are sampled from time 0 to 10 with stepsize
0.1, and the time-course curves are illustrated in
Fig.3 [10]. For this system, precisions of our obtained results are also lower than those reported
in [10], especially for the second equation. This is
because Liu et al. employed a local search procedure to refine the parameter values, while in our
method no local searching process is implemented.
Meanwhile, because derivatives of X2 fluctuate in a
small interval, our method only successfully obtain
the correct connection of the second equation with
success rate 85%.
4) Inferring Results of S4: The fourth benchmark system is a 5-D system where the parameters
are listed in Tab.V [9]. To compare with the results
reported in [9], we employ a data set of 60 samples
generated from time 0 to 0.5 with stepsize 0.0083
with initial conditions X1 (0) = 0.7, X2 (0) =
0.12,X3(0) = 0.14, X4 (0) = 0.16 and X4 (0) =
0.18, and the time-course curves are illustrated in

Fig.4. Results listed in Tab. V demonstrate that for
equations 1, 2,and 5, the best obtained parameter
values are more precise than what were reported in
[9], however, inference results for equations 3 and 4
are not satisfactory. Fig. 4 illustrates that the curves
of X3 and X4 rapidly reach equilibrium states after
short increasing procedure, and consequently, our
method, based on fitting of derivatives (slopes),
could not identify the network connections well.
B. Improving Inference Results by Multiple Data
Sets
To improve the inference results of our method,
we employ multiple data sets inferring the S-system
S4. For comparison with results reported in [9], we
generate four data set with initial conditions
1) X1 (0) = 0.7, X2 (0) = 0.12,X3 (0) = 0.14,
X4 (0) = 0.16 and X4 (0) = 0.18;
2) X1 (0) = 0.7, X2 (0) = 0.12, X3 (0) = 0.14,
X4 (0) = 0.16 and X4 (0) = 0.18;

TABLE IV: Precise and Obtained Values for Parameters of Benchmark System S3.
i
1
2
3
4
5

αi
2
(1.99,2.00)
2
(1.86,1.01)
4(3.9)
(3.93,4.00)
4
(3.97,4.00)
1
(0.98,1.00)

gi1
0
(0.00,0.00)
0.5
(0.52,0.99)
0(0)
(0.00,0.00)
0
(0.00,0.00)
0
(0.00,0.00)

gi2
0
(0.00,0.00)
0
(0.00,0.00)
0.5(0.51)
(0.51,0.50)
0
(0.00,0.00)
0
(0.00,0.00)

gi3
0
(0.00,0.00)
0
(0.00,0.00)
0(0)
(0.00,0.00)
0.8
(0.81,0.80)
0
(0.00,0.00)

gi4
0
(0.00,0.00)
0
(0.00,0.00)
0(0)
(0.00,0.00)
0
(0.00,0.00)
1
(1.00,1.00)

gi5
0
(0.00,0.00)
-1
(-1.03,-1.00)
0(0)
(0.00,0.00)
0
(0.00,0.00)
0
(0.00,0.00)

βi
2
(1.99,2.00)
4
(3.8351,4.03)
4(3.9)
(3.93,4.00)
1
(0.98,1.00)
4
(3.97,4.00)

hi1
0.5
(0.50,0.50)
0
(0.00,0.00)
0(0)
(0.00,0.00)
0
(0.00,0.00)
0
(0.00,0.00)

hi2
0
(0.00,0.00)
0.5
(0.52,0.50)
0(0)
(0.00,0.00)
0
(0.00,0.00)
0
(0.00,0.00)

hi3
0
(0.00,0.00)
0
(0.00,0.00)
0.8(0.81)
(0.81,0.80)
0
(0.00,0.00)
0
(0.00,0.00)

hi4
0
(0.00,0.00)
0
(0.00,0.00)
0(0)
(0.00,0.00)
1
(1.01,1.00)
0
(0.00,0.00)

hi5
-1
(-1.00,-1.00)
0
(0.00,0.00)
0(0)
(0.00,0.00)
0
(0.00,0.00)
0.5
(0.50,0.50)

success rate
100%
10%
85%
100%
100%

TABLE V: Precise and Obtained Values for Parameters of Benchmark System S4.
i
1
2
3
4
5

αi
5
(4.91, 4.38)
10
(9.99, 9.32)
10
(9.47, 10.88)
8
(22.37, 8.33)
10
(9.92, 9.63)

gi1
0
(0.00, 0.00)
2
(2.00, 1.79)
0
(0.00, 0.00)
0
(0.00, 0.00)
0(0)
(0.00, 0.00)

gi2
0
(0.00, 0.00)
0
(0.00, 0.00)
-1
(-1.05, -1.66)
0
(0.00, 0.00)
0
(0.00, 0.00)

gi3
1
(1.00, 1.43)
0
(0.00, 0.00)
0
(0.00, 0.00)
2
(1.89, 1.90)
0
(0.00, 0.00)

gi4
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
2
(2.06, 2.06)

gi5
-1
(-1.01, -0.90)
0
(0.00, 0.00)
0
(0.00, 0.00)
-1
(-0.54, -0.75)
0
(0.00, 0.00)

3) X1 (0) = 0.1, X2 (0) = 0.12, X3 (0) = 0.7,
X4 (0) = 0.16 and X4 (0) = 0.18;
4) X1 (0) = 0.1, X2 (0) = 0.12,X3(0) = 0.14,
X4 (0) = 0.16 and X4 (0) = 0.7,
where each data set contains 15 data samples. The
obtained results are reported in Tab. VI.
By employing four data sets generated by different initial conditions, our method can perfectly
identify the network topology of S4 with success
rate 100%, and the obtained kinetic constant values
are generally better than the results reported in [9].
However, sometimes rate constants obtained by our
method are a bit worse than the compared results,
because the rate constants are generated via the
Least Square Method (LSM) in our method.
V. O BTAINING THE GRN S ETTINGS FROM
C ANDIDATE N ON - DOMINATED S OLUTIONS
Then, we select the Pareto solutions with minimum Akaike Information Criterion (AIC) value as
the ultimate GRN inference results. The AIC for a
candidate Pareto solution is computed as


2(kgi k + khi k)
1
Ji (gi , hi ) +
,
AIC(x) = log
M
M
(10)
where M is the number of data samples [16]. Based
on the results selected from obtained Pareto sets of
20 independent runs, we compute the sensitivity Sn

βi
10
(9.89, 9.57)
10
(9.99, 10.56)
10
(9.49, 9.85)
10
(24.91, 9.89)
10
(9.89, 10.66)

hi1
2
(2.02, 1.47)
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)

hi2
0
(0.00, 0.00)
2
(2.00, 2.06)
-1
(-1.04, -1.25)
0
(0.00, 0.00)
0
(0.00, 0.00)

hi3
0
(0.00, 0.00)
0
(0.00, 0.00)
2
(1.84, 1.88)
0
(0.00, 0.00)
0
(0.00, 0.00)

hi4
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
2
(0.93, 2.24)
0
(0.00, 0.00)

hi5
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
2
(2.00, 2.14)

TABLE VII: Average Index Values of Inference
Results of S4 Obtained by Four Data Sets.
i
1
2
3
4
5

TP
3
2
3
3
2

FN
0
0
0
0.05
0

TN
6.95
8
6.6
6.1
8

FP
0.05
0
0.4
0.85
0

Sn
1
1
1
0.98
1

Sp
0.99
1
0.94
0.88
1

and specificity Sp as follows:
Sn =

TP
TP + FN

(11)

Sp =

TN
TN + FP

(12)

where TP, FN, TN and FP represent the true positive
(TP), false negative (FN), true negative (TN) and
false positive (FP) predictions of the parameters.
Average values of Sn and Sp for 20 independent
runs for S4 are listed in Tab. VII. It is shown
that the sensitivity Sn is equal to 1 except for the
fourth equation (Sn = 0.98), which means that for
all equations the true network connections can be
identified with a probability approximately equal to
1; however, sometimes the specificity Sn is less than
1, which occurs when the true network topologies
plus false positive connections can achieve better
fitting errors.

success rate
85%
100%
10%
90%
100%

TABLE VI: Improved Inference Results of Benchmark System S4 by Four Data Sets.
i
1
2
3
4
5

αi
5
(4.84, 4.38)
10
(9.99, 9.32)
10
(10.92, 10.88)
8
(8.66, 8.33)
10
(9.97, 9.63)

gi1
0
(0.00, 0.00)
2
(2.00, 1.79)
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)

gi2
0
(0.00, 0.00)
0
(0.00, 0.00)
-1
(-1.00, -1.659)
0
(0.00, 0.00)
0
(0.00, 0.00)

gi3
1
(1.01, 1.43)
0
(0.00, 0.00)
0
(0.00, 0.00)
2
(1.99, 1.90)
0
(0.00, 0.00)

gi4
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
2
(2.00, 2.06)

VI. C ONCLUSION
In this paper, we presents an evolutionary multiobjective approach to inference network of Ssystem, using the SPEM decoupling the ODE system. Then, AIC is employed to select one Pareto
solution as the final inference result. Due to the
proposed multi-objective optimization model, there
is no need to preset any parameter value before the
inference MOEA is run. So, our method could be
generally applicable to various kind of GRN model.
However, we also find that this method is sensitive
to the curve of expression level, and multiple data
sets are needed to improve its performances. Future
work will be focused on inference of GRN based
on noisy data, and how to generalize the proposed
method to infer big-scale GRNs.
ACKNOWLEDGMENT
This work is partially supported by “the Fundamental Research Funds for the Central Universities (WUT: 2014-Ia-007)” and the National Science Foundation of China under grants 61303028,
61173060 and 91230118.
R EFERENCES
[1] M. I. Davidich and S. Bornholdt, “Boolean network model
predicts cell cycle sequence of fission yeast,” PLOS One,
vol. 3, no. 2, p. e1672, 2008.
[2] N. Friedman, M. Linial, I. Nachman and D. Pe’er, “Using
Bayesian networks to analyze expression data”, J. Computat.
Biol., vol. 7, nos. 3-4, pp. 601-620, 2000.
[3] A. V. Werhli, M. Grzegorczyk and D. Husmeier, “Comparative evaluation of reverse engineering gene regulatory networks with relevance network, graphical Gaussian models
and Bayesian networks”, Bioinformatics, vol. 22, no. 20, p.
2523, 2006.
[4] K. Y. Tsai and F. S. Wang, “Evolutionary optimization
with data collocation for reverse engineering of biological
networks”, Bioinformatics, vol. 21, no. 7, p. 1180, 2005.

gi5
-1
(-1.05, -0.90)
0
(0.00, 0.00)
0
(0.00, 0.00)
-1
(-0.97, -0.75)
0
(0.00, 0.00)

βi
10
(9.88, 9.57)
10
(9.99, 10.56)
10
(10.92, 9.85)
10
(10.84, 9.89)
10
(9.97, 10.66)

hi1
2
(2.03, 1.47)
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)

hi2
0
(0.00, 0.00)
2
(2.00, 2.06)
-1
(-1.01, -1.25)
0
(0.00, 0.00)
0
(0.00, 0.00)

hi3
0
(0.00, 0.00)
0
(0.00, 0.00)
2
(1.67, 1.88)
0
(0.00, 0.00)
0
(0.00, 0.00)

hi4
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
2
(1.74, 2.24)
0
(0.00, 0.00)

[5] W. Lee and Y. Hsiao, “Inferring gene regulatory networks
using a hybrid GA-PSO approach with numerical constraints
and network decomposition”, Inform. Sci., vol. 188, pp. 8099, 2012.
[6] Y. Li, S. Jin, L. Lei, Z. Pan and X. Zou, “Deciphering
deterioration mechanisms of complex diseases based on the
construction of dynamic networks and systems analysis”,
Sci. Rep., vol. 5, e9283, 2015.
[7] X. Xiao, W. Zhang and X. Zou, “A new asynchronous
parallel algorithm for inferring large-scale gene regulatory
networks”, Plos One, vol. 10, no. 3, p. e0119294, 2015.
[8] D. Tominaga, N. Koga, and M. Okamoto, “Efficient numerical optimization algorithm based on genetic algorithm for
inverse problem,” in Proc. Genet. Evol. Computat. Conf.,
vol. 251. 2000, p. 258.
[9] L. Palafox, N. Noman, and H. Iba, “Reverse engineering
of gene regulatory networks using dissipative particle swarm
optimization“, IEEE Trans. Evol. Comp., vol. 17, no. 4, 577587.
[10] L. Liu, F. Wu, and W. Zhang , “Inference of Biological
S-system Using Separable Estimation Method and Genetic
Algorithm”, IEEE/ACM Trans. on Comput. Bio. Bioinform.,
vol. 9 no. 4: 955-965, 2012.
[11] C. Spieth, F. Streichert, N. Speer and A. Zell, “Multiobjective model optimization for inferring gene regulatory
networks”, in Proc. Evol. Mult.-Criter. Opt. Conf., vol. 3410,
2005, pp. 607-620.
[12] P. Liu and F. Wang, “Inference of biochemical network
models in s-sytem using multi-objective optimization approach”, Bioinformatics, vol. 31, no. 23, pp. 1085-1092,
2006.
[13] P. Koduru, Z. Dong, S. Das, S. M. Welch, J. L. Roe and
E. Charbit, “A multiobjective evolutionary-simplex hybrid
approach for the optimization of differential equation models
of gene networks”, IEEE Trans. Evol. Comp., vol. 12, no.
3, pp. 572-590, 2008.
[14] X. Cai, Z. Hu, S. Das and S. M. Welch, “A hierarchical
Pareto dominance based multi-objective approach for the
optimization of gene regulatory network models”, in Proc.
Evol. Comput. Conf., 2012, pp. 1-6.
[15] Y. Chen, W. Xie and X. Zou, “A binary differetial evolution
algorithm learning from explored solutions”, Neurocom.,
vol. 149, 1038-1047, 2015.
[16] C. .W. Yao B. D. Hsu and B. .S. Chen, “Constructing
gene regulatory networks for long term photosynthetic light
acclimation in Arabidopsis thaliana”, BMC Bioinform., vol.
12, e335, 2011.

hi5
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
0
(0.00, 0.00)
2
(2.01, 2.14)

success rate
100%
100%
100%
100%
100%

