Fuzzy-Based Self-Interactive Multiobjective Evolution Optimization for Reverse Engineering of Biological Networks

Fuzzy-Based Self-Interactive Multiobjective

Evolution Optimization for Reverse Engineering

of Biological Networks

Shinq-Jen Wu, Cheng-Tao Wu, and Jyh-Yeong Chang, Member, IEEE

Abstract—S-system modeling from time series datasets can

pro-vide us with an interactive network. However, system identifica-tion is difficult since an S-system is described as highly nonlin-ear differential equations. Much resnonlin-earch adopts various evolution computation technologies to identify system parameters, and some further achieve skeletal-network structure identification. However, the truncated redundant kinetic orders are not small enough as compared with the preserved terms. In this paper, we integrate quantitative genetics, bacterium movement, and fuzzy set theory into evolution computation to develop a new genetic algorithm to achieve convergence enhancement and diversity preservation. The proposed exploration and exploitation genetic algorithm (EEGA) can improve the best-so-far individual and ensure global optimal search at the same time. The EEGA enhances evolution conver-gence by golden section seed selection, normal-distribution repro-duction, mixed inbreeding and backcrossing, competition elitism, and acceleration operations. Search-then-conquer evolution direc-tion operadirec-tions, eugenics-based screen-sifting mutadirec-tion, eugenic self-mutation, and fuzzy-based tumble migration preserve popu-lation diversity to avoid premature convergence. Furthermore, to ensure that a reasonable gene regulation network is inferred, fuzzy composition is introduced to derive a reconstruction index. This performance index let EEGA possess self-interactive multiobjec-tive learning. The proposed fuzzy-reconstruction-based multiob-jective genetic algorithm is examined by three dry-lab biological systems. Simulation results show that a safety pruning action is guaranteed (the truncation threshold is set to be 10−15), and only one- or two-step pruning action is taken.

Index Terms—Multiobjective, real-value coding, self-interactive,

structure identification.

I. INTRODUCTION

T

re-Manuscript received March 21, 2011; revised August 12, 2011 and December 8, 2011; accepted December 19, 2011. Date of publication February 7, 2012; date of current version October 2, 2012. This work was supported by the National Science Council, Taiwan, under Grant NSC-99-2221-E-212-021.

S.-J. Wu is with the Department of Electrical Engineering, Da-Yeh University, Changhua 51591, Taiwan (e-mail: jen@mail.dyu.edu.tw).

C.-T. Wu and J.-Y. Chang are with the Department of Electrical and Control Engineering, National Chiao-Tung University, Hsinchu 300, Taiwan (e-mail: dau@cn.nctu.edu.tw; jychang@mail.nctu.edu.tw).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TFUZZ.2012.2187212

peated modification and an undue amount of experiment data are necessary in parameter identification, especially for a system with many substances or reactions involved. S-system struc-ture [3], [4] is another preferred nonlinear dynamic model. This model can uniquely map dynamic interaction onto its parameters and possesses good generalization characteristics. Some researchers used gradient-based computation technolo-gies to identify the parameters of an S-system model. Marino and Voit [5] wrote an algorithm to gradually increase model complexity. Chou et al. [6] adopted an alternating regression (AR) method. Vilela et al. [1] further proposed an eigenvector optimization method to solve the convergence issues of an AR approach. Kutalik et al. [7] introduced the Newton flow analysis. Recently, many researchers have tried to infer a gene regula-tory network by stochastic search intelligent technologies such as genetic programming [8]–[11], evolutionary algorithms [12], evolution strategies [13], differential evolution [14]–[18], ge-netic algorithms (GAs) [19]–[21], simulated annealing [22], hybrid GA and simulated annealing [23], radial basis function networks [24], and a neural network with particle swarm opti-mization [25], [26]. However, as system states and parameters increase, parameter estimation becomes difficult and the solu-tions are problematic. Furthermore, pruning redundant kinetic orders to infer a suitable network structure is a big challenge. Based on the fact that a genetic network is connected sparsely [27], various penalty terms are introduced [12], [14], [19], [28]. However, no matter what kinds of penalty terms are chosen, the truncated redundant kinetic orders are not small enough to guarantee that a safely pruning operation is taken.

Technological contributions of this paper can be described as follows. A new GA with various advanced genetic op-erators is proposed in Section II-A. The proposed GA pos-sesses a tradeoff between exploration and exploitation. In Section II-B, single-objective performance candidates for pa-rameter identification are chosen from 28 indexes. These can-didates are fuzzily integrated with a new kinetic order index for multiobjective structure identification. In Section III, the proposed fuzzy-reconstruction-based multiobjective genetic al-gorithm (FRMOGA) is examined by three genetic networks. Section IV is the conclusive remark.

II. FUZZY-BASEDSELF-INTERACTIVEMULTIOBJECTIVE EVOLUTIONCOMPUTATION

S-system is a well-known canonical nonlinear model to cap-ture genetic interactions/transcriptional regulation. Based on 1063-6706/$31.00 © 2012 IEEE

Fig. 1. Golden section.

biochemical system theory, the net influx (V_i+) and efflux (V_i−) of a system are approximated as power law functions. Each individual metabolite, protein, or gene is described as

˙ Xi= Vi+− Vi−= αi n + m j = 1 Xgi j j − βi n + m j = 1 Xhi j j for i = 1, 2, . . . , n (1) where n and m are the numbers of dependent and independent variables, respectively; αi and βj are rate constants, gij and

hij are kinetic orders to denote the interaction from Xj to Xi,

where a positive value denotes excitatory effect and negative for inhibitory effect. Recently, various evolutionary optimiza-tion technologies were used to infer S-type gene regulatory networks (parameter estimation and structure identification of an S-system). However, as system states increase, the inference of a suitable network structure becomes a big challenge and the solutions are problematic. Local minimum solutions will generate an error structure. Avoiding from sticking into local minima is very important to infer such a high dimensional and nonlinear system by computation approach. Sometimes, in a bi-ological system, an interaction order is small, but the interaction effect cannot be neglected. Besides, noise and uncertainty exist in a biological system, especially for dealing with microarray data. To achieve correct structure, the gap between true and redundant interactions should be increased. In other words, the pruning threshold should be small enough to ensure safety prun-ing. Furthermore, a wide searching space is necessary since no prior information is provided. Therefore, the used technologies should possess good search power and has the ability to escape from local minima. It is hard for a state-of-the-art GA with sim-plex crossover (SPXGA) to obtain a satisfactory solution in a limited computation time [21]. Ho et al. proposed a GA with in-telligent crossover and Cauchy–Lorentz probability-distribution mutation to identify genetic networks [21]. We here propose an exploration and exploitation genetic algorithm (EEGA) to improve the best-so-far individual and to ensure global optimal search at the same time. EEGA enhances evolution convergence by golden section seeds selection, normal-distribution reproduc-tion, mixed inbreeding and backcrossing, competition elitism, and acceleration operations. Search-then-conquer evolution di-rection (SCED) operation, eugenics-based screen-sifting mu-tation, eugenic self-mumu-tation, and fuzzy-based tumble

migra-so-far individual to further improve its local and global search power.

1) Population Initialization: We adopt real coding to exploit the gradualness of continuous variables. Each parameter is ini-tialized to cover an entire search space. Therefore, an initial population with NP individuals is chosen as

Ii0 = Im in+ r(Im ax− Im in), i = 1, ..., NP (2)

where r is a real number from an uniform distribution [0,1], Im in

and Im ax are the lower and upper bounds of individual Ii.

2) Golden Section Seeds Selection: Individuals in a popula-tion are arranged according to their fitness values in a descending order. Then, a golden section method in Fig. 1 is used as follows to choose three Elitism-seeds groups, i.e., n1S, n2S, n3S:

n1S = [τ· nS] for the first seeds group

n2S = [τ· (nS − n1S)] for the second seeds group

n3S = nS − n1S − n2S for the third seeds group (3)

where [•] is a Gauss mark, and τ is the golden section con-stant; the number of selected seeds is nS = [τ·NP] for a small

population (NP < 10) and nS = [(1− τ)·NP] for NP ≥ 10.

3) Search-Then-Conquer Evolution Direction: We then choose three preferable individuals (Ib, Is, and It) from these

three seeds groups; one individual from each group. Fb, Fs,

and Ft are the associated fitness. These three individuals are

used in SCED operation to determine the evolutionary direction of a population. SCED proceeds a three-phase search. In the early-evolution phase, a random walk search is taken for wider searching

Isced = Ib+ r1∗ D1∗ (Ib− Ir1) + r2∗ D2∗ (Ib − Ir2). (4) In the middle-evolution phase, a directed-walk search is taken to provide a good searching direction

Isced = Ib + r1∗ D1∗ (Ib− Is) + r2∗ D2∗ (Ib− It). (5)

Then, a conquer strategy is used to achieve fast convergence for the final phase

Isced = Ib+ r1∗ D1∗ (Ib− Ib1) + r2∗ D2∗ (Ib− Ib2) (6)

where r1 and r2 are random factors; D1 and D2 denote

evo-lutionary directions; Ir1 and Ir2 are two randomly selected individuals; and Ib1 and Ib2 are the perturbed values of Ib.

The associated fitness value Fsced is estimated. If Fsced is

better than one of the three preferable fitness values, the in-dividual is replaced. Fig. 2 describes the SCED operation. SCED has the ability for local and global search synchronously. We can use this operator to achieve widely searching and

Fig. 2. SCED.

then quickly converging toward a global optimal solution. The operator helps us from bogging down into a local optimal solution.

4) Normal-Distribution Reproduction: Let the fitness of off-spring be a Gaussian distribution. [0.5σNp] individuals are

re-lated to the first seeds group only, [σNp] individuals are related

to the first two seeds groups only, and [1.5σNp] individuals

are related to these three seeds groups only. For diversity, the residue individuals are chosen from those eliminated individ-uals only. σ is the standard deviation of a normal distribution. Fig. 3 describes a parent distribution for crossover. During a crossover operation, if the same individual is chosen for a father and a mother, the individual is perturbed with a Gaussian/normal function N (0, σ). Therefore, we denote the groups in Fig. 3 by

Fig. 3. Normal-distribution reproduction.

perturbed groups, where Iisdenotes the individuals from the ith

seeds group ˜

Fig. 4. Mixed inbreeding and backcrossing operation.

5) Competition Elitism: Competition carries on in each ge-netic operation. Soft competition is used for evolution direction operation and mutation. Three seeds groups, i.e., I1s, I2s, I3s,

take over the most responsibility; the more highly fit strings have a higher number of offspring. Only the mutated genes with better fitness values than their precursors can be accepted. Winner-take-all is adopted for acceleration, self-mutation, and migration. Only the best gene is chosen to participate in these operations. Competition also occurs between generations; only the winner can survive to the succeeding generation.

6) Mixed Inbreeding and Backcrossing: Randomly choose two points in parent strings (A and B) for two-point crossover to generate children C and D. These two children or just one is selected to replace the father or the mother for further crossover. In other words, the selected candidates are two children (C and D) for close inbreeding, or a parent and a child (A and D or B and C) for backcrossing. Then, a two-point crossover operation is used again to generate offspring E and F. Fig. 4 describes the proposed mixed inbreeding and backcrossing operation.

The inbreeding coefficient of offspring X1from two children

(C and D) is denoted by fx1, and that of offspring X2 from backcrossing between parent and children (A and D or B and C) is denoted by fx2 [29] fx1 = fC D = 1 4(fA A+ 2fA B + fB B) = 1 4  1 +1 2(FA+ FB) + 2fA B  fx2 = fA D = 1 2(fA A+ fA B) = 1 2  1 2(1 + FA) + fA B  (8) where FAand FB are the inbreeding coefficient of individuals,

A and B, respectively. The average relative coefficient fav with

backcrossing probability rb is defined as

fav = rb ·

1

4(fA D+ fA C + fB C + fB D) + (1− rb)· fC D. (9) Therefore, fav= fav,m axfor parent coming from the same seeds

group and fav = fav,m infor no blood relationship parent. Based

on the average relative coefficient, we set the selected candidates

Fig. 5. Eugenics-based screen-sifting mutation.

for further crossover to be B and C in case of rb≤ fav,m in, to be

A and D for fav,m in< rb ≤ fav,m ax, and to be C and D for rb >

fav,m ax.

7) Eugenics-Based ScreSifting Mutation: In order to en-sure global search and improve convergence, we here propose a mutation-then-eugenics operation. In stead of adopting only one gene (one-point mutation), two genes (two-point mutation), or some fixed-genes (mask mutation) for mutation, we let all genes or randomly choose some genes to mutate with mutation probability assigned by a designer. Those genes with qualified mutation rate will mutate (screen sifting). After this kind of mu-tation, a variety population is generated. We wish the mutated population distributes over a wide but not an entire space. A deviation threshold is then defined as

T A =|Fmut− Fb| Fb

(10) where Fmutis the fitness value of the mutated individual. Then,

the mutated individuals compete with their source individuals. Those winners or losers with acceptable deviation (TA < 10) can pass down. In engineering, a ratio 3 or 5 means much far away or apart in space. Here, we use 10 to denote that mutated individuals distribute widely. The proposed eugenics-based mu-tation in Fig. 5 can not only increase population diversity but can guarantee fine offspring as well.

8) Acceleration: A steepest descent method or simplex method is used to enhance convergence. After acceleration, the best individual is Ib =  ¯_{Ib, if f ( ¯Ib) < f (I}0 b) ¯ Ib− α∇f, otherwise (11) where I0

b and ¯Ib are the best individuals before and after both

mutation and crossover operations, α is the size of a step to determine the descent rate, and∇f is the gradient of an objective

function. In other words, an acceleration operation starts when both crossover and mutation operations can no longer increase fitness.

9) Eugenic Self-Mutation: Inherent premature convergence and leak mountain climbing are two weak points of real-coding GA. Acceleration operations can compensate for the latter. We here propose eugenic self-mutation to avoid premature con-vergence. The mutation is applied to the best individual only. Each gene in the best individual is chosen and assigned a mu-tation probability. Competition carries on between the original

individual and the mutated individual. The winner will pass down to the next operation.

10) Fuzzy-Based Tumble Migration: We here integrate fuzzy concept and the movement of an E.coli bacterium into a migration operation. A population P is composed of Np

in-dividuals (I1, I2,. . ., INP). An individual Ii is composed of n

chromosomes (xi1, xi2, . . ., xin): P ={Ii}i= 1,...,NP = ⎡ ⎢ ⎢ ⎣ I1 I2 .. . INP ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ x11 x12 . . . x1n x21 x22 . . . x2n .. . xNP,1 xNP,2 . . . xNP,n ⎤ ⎥ ⎥ ⎦ Ii={xij}_{j = 1,...,n}= [ xi1 xi2 . . . xin] . (12)

Define fuzzy sets ˜Ii, i= 1, . . . , n, to represent individuals not

close to Ib. The relative cardinality of ˜Iidescribes the degree of

the deviation of Iifrom Ib:

|˜Ii|rel= _|U|˜Ii|

i| ˜ Ii=  (xij, μI˜i(xij))∀j = 1, . . . , n = [ μi1 μi2 . . . μin] μij= μ_I˜_i(xij) =  xij− xbj xbj   (13)

where Ui are the universal set of Ii. Therefore, the following

fuzzy set ˜P denotes the deviation tendency of a population to the best individual Ib, and its α level cut ˜Pα contains the elements

of a universal set U with membership grades greater than or equal to α (closeness tolerance). In other words, ˜Pα contains

the elements in ˜P far enough to Ib. Therefore, we define the

population diversity of P as η: η=Δ   ˜Pα |U| ˜ P = ⎡ ⎢ ⎢ ⎢ ⎣ ˜ I1 ˜ I2 .. . ˜ IN P ⎤ ⎥ ⎥ ⎥ ⎦= ⎡ ⎢ ⎢ ⎣ μ11 μ12 . . . μ1n μ21 μ22 . . . μ2n .. . μNP,1 μNP,2 . . . μNP,n ⎤ ⎥ ⎥ ⎦ . (14) If the population diversity η is smaller than a diversity threshold

ξ, then migration starts. The tumble movement of an E.Coli cell

is adopted for migration operation: After a tumble migration, a new individual is generated as

Inew = Ib+ c·  Δ √ ΔT_Δ  (15) where c is the size of the step taken in a random direction, and Δ indicates a vector in a random direction.

B. Preference-Based Self-Interactive Multiobjective Optimization

To infer a physically realizable skeletal-network structure, various pruning-penalty performance indexes are proposed. On the basis of that a genetic network is connected sparsely [27], Kikuchi et al. [19] introduced kinetic orders as the pruning (penalty) terms of a fitness function. Kimura et al. [12] replaced kinetic orders by their descending orders and introduced a maxi-mum in degree I to determine the maximaxi-mum number of genes that affect the ith gene. Norman and Iba [14] modified penalty terms by the sorted orders of both excitatory and inhibitory terms. Liu and Wang [28] introduced ε-constrain weighted sum. Ko et al. [30] considered a flux connectivity relationship and in-troduced a sensitivity term. However, no matter what kinds of penalty terms are used, pruning approaches face a challenge: Weighting coefficients need to be carefully tuned, and there are no clear guidelines for setting suitable penalty weights [31]. Furthermore, even if so much effort has been done, the pruned criteria in Table V, shown in the Appendix, are not small enough as compared with the preserved terms. In other words, a safely pruning action is not guaranteed. In the previous section, we integrate quantitative genetics, bacterium movement, and fuzzy set theory into evolution computation to develop a new GA to achieve convergence enhancement and diversity preservation. However, a performance index that is defined on various errors is a criterion to show how closing the estimated data to a real profile. This index determines both search directions and com-putation time. Therefore, the choice of a performance index is a key point for computation optimization. In this section, we shall get single-objective performance index candidates first. Those candidates are further integrated into a self-interactive multiobjective function.

1) Single-Objective Performance Index Candidates: The proposed EEGA algorithm with 28 performance indexes in Table VI, shown in the Appendix, is used to identify the pa-rameters of S-systems for two dry-lab experiments in Figs. 6(a) and 6(b). Those well-performance indexes are chosen as candi-dates, and are further integrated into a multiobjective function. All computation is performed on an Intel core duo 3.16-GHz computer using Microsoft Windows XP. Parameter ranges are [0, 30] for rate constants and [–4, 4] for kinetic orders. 250 000 maximum iterations are performed to minimize those error cri-teria. Eight sets of experiment data xiexp and the corresponding

slope information ˙xi

expare generated from the dynamic system

in [20] for a genetic branch pathway, and in [15] for a cascade pathway. The simulation time of an experiment is from t = 0 s to t = 8 s, and sample time is set to be 0.02.

The genetic branch pathway in Fig. 6(a) has four dependent constituents, i.e., x1, x2, x3, and x4, and one constant source,

i.e., x0. The results in Table VII, shown in the Appendix, give us

some information. First, summation operations perform much better than the corresponding maximum operations except for

Jl

5, J7l, and J14l . Second, J1l, J2l, and J3l are, respectively, better

than Jl

4, J5l, and J6l; J7l is better than J13l ; and J10l is

bet-ter than Jl

14. Those results show that even if errors are

Fig. 6. Three gene regulation systems and their S-type Models. (a) Genetic branch pathway with two regulatory signals [20]. (b) Cascade pathway with three steps and two feedback signals [15]. (c) Small-scale genetic network [32], [19].

Fig. 7. Block diagram of an artificial system for structure identification. down errors to generate a small-quantity criterion for updat-ing, this operation still improves performance. Third, indexes with normalization have better performance (Jl

10, J11l , and J12l

are, respectively, better than Jl

7, J8l, and J9l). Fourth, among

summation operations, the performance of those weighting-related indexes is always ta > 1 > ts (“>“ denotes better).

However, if slope normalization is taken (Jl

10, J11l , J12l ), then

it becomes 1 > ta > ts. Based on these four phenomena, we

suggest to choose JS

1, J2S, J3S, J9S, J10S, J12S , J3M, and J10M as

learning-criterion candidates.

We now consider another biological system. The cascade pathway in Fig. 6(b) has three dependent constituents x1, x2,

and x3, and one constant sourcex4. Table VIII in the Appendix

has similar results as in Table VII except the third phenomenon. This is due to a dramatically sharp-slope point in x2 profile.

Therefore, slope normalization is not suitable. Based on Ta-ble VIII, shown in the Appendix, we know that J1S, J3S, J6S, J7S,

JS

9 , J13S , J1M, J3M, and J9M are all suitable for learning indexes.

In the following section, we shall use JS

3 and J10S as our

error-and slope-performance cerror-andidates to develop a multiobjective optimization technology.

2) Fuzzy-Reconstruction-Based Multiobjective Genetic Al-gorithm: In the previous section, we examine parameter

iden-tification performance under 28 criteria. The results give us a very important concept. Except normalization operation, the dynamic behavior of a system is an important factor. We now extent the results to define a realizable learning index for the skeletal-structure identification of a biological system. Concen-tration error Je and kinetic order Jo are two key penalties to

examine fitness between measured data and estimated data, and to ensure that a genetic network is sparsely connected. Besides these two, we also introduce a slope error penalty Jd to achieve

smooth time evolution. Furthermore, at any time instant, the net interaction strength (row sum of matrix K) of all species to some species is finite. Therefore, we define a kinetic order penalty as

J0 =K_∞= max i ⎧ ⎨ ⎩ n  j = 1 |kij| ⎫ ⎬ ⎭ (row sum) K ={kij, i = 1, ..., n, j = 1, . . . , 2(n + m)} = ⎡ ⎢ ⎣g11 . . . g1,n + m .. . h11 . . . h1,n + m .. . ... gn 1 . . . gn ,n + m ... hn 1 . . . hn ,n + m ⎤ ⎥ ⎦ . (16) JS

3 and J10S are chosen as our error- and slope-penalty,

respectively.

As we know, concentration error Je, slope error Jd, and

skeletal-structure penalty Jo are belong to different scales. We

here normalize them as Je = Je/Jee, Jd = Jd/Jde, and J0 =

Jo/Joe, where Jee, Jde, and Joe are the expected concentration

error, slope error, and interaction measure, respectively. To infer a biological regulation network, our objective is to push both Je

and Jd to approach zero but to obtain a nonzero minimum value

of Jo. The targets of Jeand Jdare the same, but they are totally

different from that of Jo. Therefore, summating their scaled

functions is not suitable. Moreover, the structure identification of a biological system is to get minimum value of Jo under

TABLE I

VARIABLES ANDMODULUS OFSUBPROCESSES INFIG. 7

allowable Jd and Je. It is different from fuzzy rule selection,

where rules are extracted only if it is possible to either maintain or even improve the system’s accuracy [33], [34] and general multiobjective optimization problems. Therefore, Pareto-optimal-based evolution computation approaches [35]–[41] are also unsuitable. In this paper, we use fuzzy composition “” (min max) to integrate these three penalties into a reconstruction performance. Based on this performance, the proposed GA becomes a self-interactive multiobjective evolution optimization technology.

Biological systems are dynamic changeable systems. To guar-antee well performance of computation approaches for gene net-work identification, except concentration error, we should con-sider slope error for data with close values but different slope. Additionally, sparse connection and the correct division of exci-tory and inhibiexci-tory interactions are two important issues that are to be concerned. Therefore, the whole structure identification can be physically realized as three modules in series (a state estimator, S-system modeling, and a slope estimator). Since un-certainty and noise are two serious issues in a biological system, fuzzy relation and fuzzy functions are adopted to describe the dynamic behavior of subprocess modules. Fig. 7 shows a block diagram to describe the inference of structure identification.

Rx, R, and Rd are, respectively, fuzzy modules for a state

estimator, an S-system, and a slope estimator. R−1is the inverse fuzzy relation for S-system modeling. State estimator Rx(x, ˆx)

describes the fuzzy relation between input x (true state) and output ˆx (estimated state). The corresponding fuzzy relation

equation is ˆx = xR, where “” denotes min max composi-tion. Here, our purpose is to minimize Je such that estimated

state ˆx approaches to true state x. Similar operations are used

for S-system modeling R(ˆx, ˙ˆx) and slope estimator Rd( ˙ˆx, ˙x).

Table I lists the inputs, outputs, fuzzy equations, and optimiza-tion of these subprocesses. Steps 1–4 describe the inference flow for structure identification. Three normalized errors, i.e.,

Je, Jd, and J0, denote the distortion of these subprocess. The

entire process is composed of a forward process (RxR) and

a backward process (RdR−1). Therefore, we have

x = ˙xR−1(backward S-system modeling) = ( ˙ˆxRd)R−1(slope estimator)

= ((ˆxR)Rd)R−1(forward S-system modeling)

= (((xRx)R)Rd)R−1(state estimator)

= xRxRRdR−1. (17)

For perfect construction, the membership functions of compos-ite fuzzy relation R∗= Rx  R Rd  R−1always belong to 1.

Therefore, the reconstruction ability of structure identification is dependent on the composition of the distortion of subsystems Je, Jo, and Jd. Therefore, we define a reconstruction

perfor-mance Jrec−0 as the union (U) of JeJ0 and JdJ0 (forward

AND backward processing) min

k Jrec-o = mink [(JeJ0)U(JdJ0)]

= min k [(JeU Jd)J0] = min k [mink (Je∨ Jd)∨ J0] = min k [Je∨ Jd ∨ J0] (18)

where maximum∨ is adopted as the union operation. In other words, Jrec-0 = Je∨ Jd∨ J0.

Based on this definition, an automatic self-interactive learning in Table IX, shown in the Appendix, is done. It is noticed that Jee and Jde are values approaching zero, but Joe is a value

between 0.001 and 10. As J0 is comparable with Je or Jd,

Je and Jd are values approaching zero. In the beginning of a

learning process, concentration error Jeand slope error Jd are

much larger than J0. During the mid-learning stage, those three

penalties compete with each other to get the opportunity to be reduced. At the final stage, Jo becomes the only winner, since

both Je and Jd are close to zero. Therefore, choosing Jrec-0

can divide our learning into two phases. In phase I, the learning purpose is to get an allowable space; concentration error Je is

the key learning index and slop error Jdis for smooth evolution.

In phase II, J0is the only kernel index to achieve spare network

structure. Therefore, we modify Jrec-0 by adding phase I and

phase II weighting factors wI and wI I, i.e., Jrec-1 = max{wI I

Je, wIJd, J0}. This modification makes the performance index

suitable for all kinds of biological systems (wI = wI I = 1

for most systems). Besides, the dynamic behavior of a system is also an important factor for time-series-data-based learning. (We get this result from parameter identification in Section II-B1.) Considering dynamic behavior helps us to simultaneously fine-tune the structure and parameters of a system in phase-II learning. Therefore, we further modify the performance as Jrec= max{ηI I Je, wI Jd, J0} with an adaptive-learning ratio

ηI I = wI I × μ, where μ is an adaptive dynamic factor.

The reconstruction index is embedded into EEGA to infer biological regulatory networks from time series data. Fig. 8

Fig. 8. Logical flow of FRMOGA for skeletal-network structure identification.

describes the proposed FRMOGA for skeletal-network structure identification.

III. S-SYSTEMMODELING OFGENEREGULATIONNETWORKS We now examine the proposed FRMOGA by identifying three biological systems (genetic branch pathway, cascade pathway, and small-scale genetic networks). Search space is set to be [0, 30] for rate constants and [–4, 4] for kinetic orders. In this paper, no further assumption for self-interaction (gii and hii)

is made. All dependent and independent variables are included in an S-system. In other words, our learning is based on the super structure of an S-system model, where there are 2n(n + m+ 1) parameters to be identified. Table V in the Appendix lists the assumptions, dataset, and pruning thresholds used in the published papers. We use a cubic spline technology to get a smooth profile for time series data. Then, an integral-based-modified collocation method with piecewise linear Lagrange

polynomials as shape polynomials is adopted to approximate the generated dynamic profile [18]

xi(tl) = xi(tl−1) + 0.5ηl{fi(xi(tl), θ) + fi(xi(tl−1), θ)}

i = 1, ..., n, l = 1, ..., NS (19)

where xi(tl) and fi(xi(tl), θ) are the expansion coefficients of

the ith state and rate functions at the lth collocation point, and ηl

is a time interval. Decoupled technology is adopted to get initial values for entire coupled S-type system learning.

A. Genetic Branch Pathway System

Eight sets of concentration data are generated from a true S-system; parameters are listed in Step 0 in Table II. Simulation time is set to be 8 s and sample time is 0.02. In this system, there are 48 parameters that are to be identified. Table II shows

TABLE II

TRUE ANDESTIMATEDPARAMETERS OF ANS-TYPESYSTEM FOR ABRANCHPATHWAYNETWORK INFIG. 6(a)

TABLE III

Fig. 9. Convergence results of EEGA, HDE [15], [17], [18], improved GA [42], and DE are examined by a wide search space ([0, 100] for rate constants and [–100, 100] for kinetic orders) with a bad initial start (80 for all parameters). Solid curves represent the estimated profiles of EEGA and HDE; dashed curves represent those of the improved GA (GA+_{) and DE. Those curves are drawn for fitness evaluation from 1000 to 50 000.}

that only two-step pruning action is taken. An obvious value gap exists between redundant interactions and possible interac-tions in Step 1. We subtract redundant interacinterac-tions denoted by underlines with a threshold 10−15. The pruned structure is fur-ther learned by FRMOGA to get the modified structure in Step 2. Step 3 shows the finalized structure. The inferred structure

is identical to true structure (Step 0), and parameter values are nearly the same as those in the true system.

B. Cascade Pathway System

Eight sets of concentration data are generated from true equa-tions with parameters in Step 0 in Table III. Simulation is from

Fig. 10. Convergence comparison of EEGA with IGA and SPXGA cited from [21, Fig. 2] for a small-scale (five genes), a medium-scale (20 genes), and large-scale (30 genes) genetic networks. Six-set time series data with 11 sample points are used.

time t = 0 s to t = 8 s with a sample time of 0.02. The cor-responding super structure of the S-system has 30 parameters that are to be identified. Step 1 in Table III shows the first-time learning result. Those interactions with kinetic order less than 10−15 are deleted. The pruned model is further modified into Step 2 in Table III by FRMOGA. After retruncating trivial terms and performing learning again, we have the final structure and associated parameters in Step 3 in Table III.

C. Small-Scale Genetic Network

Our third case is a small-scale genetic network with two regulatory signals. True rate constants and kinetic orders are listed in Step 0 of Table V in the Appendix. Each experiment is simulated from time t = 0 s to t = 0.5 s with a sample time of 0.0125. In this system, there are 90 parameters that are to

be identified. Inferred parameters and kinetic orders are shown in Steps 1 and 2. Pruning threshold is also set to be 10−15. Those redundant connections to be truncated are denoted by underlines. In this system, only two-step learning is taken.

D. Convergence Comparison

Figs. 11–13, shown in the Appendix, describe the conver-gence of each step of EEGA for those three systems. Table V gives the comparison of the used dataset, assumptions, S-system models, and pruning thresholds for EEGA and the published computation algorithms. Pruning ratio ρ is defined as the ra-tio of the smallest preserved term to pruning threshold. A safety pruning is taken by EEGA with ρ = 7.5× 1012_{. Except mutation}

to escape from local minima. Some advanced operations on the best-so-far individual are to accelerate convergence. To show the global search power, a bad initial start in a very wide range is used in Fig. 9. Fig. 9 shows the convergence comparison of EEGA with DE, HDE [15], [17], [18] and our previously work improved GA (GA+_{) [42] for low- and medium-dimensional}

systems (N = 3, 4, and 20 genes) by a wide search space ([0, 100] for rate constants and [–100, 100] for kinetic orders) with a bad initial start (80 for all parameters). Fig. 10 shows con-vergence results of EEGA, intelligent two-stage evolutionary algorithm (IGA) [21], and SPXGA [19] for low-, medium-, and high-dimensional systems (N = 5, 20, and 30 genes) with the same search region as [21] ([0, 15] for rate constants and [–3, 3] for kinetic orders).

IV. CONCLUSION

The inverse problem of identifying a dynamic biological sys-tem from time series data is a central theme in syssys-tems biol-ogy. However, even if multiobjective computation approaches are used, structure identification is still a big challenge. How to avoid from sticking into local minima is very important to infer such a high dimensional and nonlinear system by com-putation approach. Additionally, a big gap between true and redundant interactions is necessary to ensure a correct

struc-ture, since noise and uncertainty exist in a biological system. In other words, the pruning threshold should be small enough. The used technology should possess good search power and has the ability to escape from local minima. In this paper, we propose an EEGA learning technology to achieve the improve-ment of the best-so-far individuals and to ensure global optimal search. A fuzzy-based multiobjective construction performance index is derived to infer a physically reasonable regulation net-work. Simulation results in Tables II–IV show that a very big gap exists in true and redundant interactions. The truncation threshold can be decreased down to 10−15, and a correct struc-ture is achieved after only one- or two-step pruning operations. Table V in the Appendix shows our pruning action is absolutely safe as compared with that in the published papers. Fig. 9 shows the good global search power of EEGA with a bad initial start (80 for all parameters) in a wide search space ([0, 100] for rate constants and [–100, 100] for kinetic orders) for three-, four-, and 20- dimensional systems. In this paper, we did not add noise in our time series data. Most researches use existing smoothing or filtering technique to deal with noise. However, those ap-proaches depend too much on smoothing skills and lack clear guidelines, especially for structure identification. In the future, we shall propose a new estimator (filter) to solve those data se-riously contaminated by white noise and color noise, or systems with bias and uncertainty.

APPENDIX TABLE V

STRUCTUREIDENTIFICATIONCOMPARISON

TABLE VI

Fig. 11. Convergence of EEGA for each step in the structure identification of a branch network in Fig. 6(a) starting from (βi, hi i) = (30, 4). Dash line denotes true value. Black solid curve shows the convergence of Step 1, blue for Step 2, and Green for Step 3.

TABLE VII

PERFORMANCECOMPARISON FOR AGENETICBRANCHPATHWAY INFIG. 6(a)

TABLE VIII

TABLE IX

UpdatedTERMS INSELF-INTERACTIVELEARNING

Fig. 12. Convergence of EEGA for each step in the structure identification of a cascade network in Fig. 6(b) starting from (βi, hi i) = (30, 4). Dash line denotes true value. Black solid curve shows the convergence of Step 1, blue for Step 2, and Green for Step 3.

Fig. 13. Convergence of EEGA for each step in the structure identification of a small-scale genetic network in Fig. 6(c) starting from (βi, hi i) = (30, 4). Dash line denotes true value. Black solid curve shows the convergence of Step 1 and blue for Step 2.

ACKNOWLEDGMENT

The authors would like to thank Prof. F.-S. Wang of the Chem-ical Engineering Department, National Chung-Chen University, for his help with reverse engineering technologies.

REFERENCES

[1] M. Vilela, I. C. Chou, S. Vinga, A. T. R. Vasconcelos, E. O. Voit, and J. S. Almeida, “Parameter optimization in S-system models,” BMC Syst. Biol., vol. 2, no. 35, 2008.

[2] L. Michaelis and M. I. Menton, “Die kinetik der invertinwirkung,” Biochem. Z., vol. 49, pp. 333–369, 1913.

[3] M. A. Savageau, Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology. Reading, MA: Addison-Wesley, 1976.

[4] E. O. Voit, Computational Analysis of Biochemical Systems: A Practical Guide for Biochemists and Molecular Biologists. Cambridge, U.K.: Cambridge Univ. Press, 2000.

[5] S. Marino and E. O. Voit, “An automated procedure for the extraction of metabolic network information from time series data,” Bioinform. Comput. Biol., vol. 4, no. 3, pp. 665–691, 2006.

[6] I. C. Chou, H. Martens, and E. O. Voit, “Parameter estimation in bio-chemical systems models with alternating regression,” Theor. Biol. Med. Model, vol. 3, no. 25, 2006.

[7] Z. Kutalik, W. Tucker, and V. Moulton, “S-system parameter estimation for noisy metabolic profiles using Newton flow analysis,” IET Syst. Biol., vol. 1, pp. 174–180, 2007.

[8] E. Sakamoto and H. Iba, “Inferring a system of differential equations for a gene regulatory network by using genetic programming,” in Proc. Congr. Evol. Comput., 2001, vol. 1, pp. 720–726.

[9] S. Ando, E. Sakamoto, and H. Iba, “Evolutionary modeling inference of gene network,” Inf. Sci., vol. 145, pp. 237–259, 2002.

[10] D. Y. Cho, K. H. Cho, and B. T. Zhang, “Identification of biochemical networks by S-tree based genetic programming,” Bioinform., vol. 22, pp. 1631–1640, 2006.

[11] H. Wang, L. Qian, and E. Dougherty, “Inference of gene regulatory net-works using S-system: A unified approach,” IET Syst. Biol., vol. 4, pp. 145–156, 2010.

[12] S. Kimura, K. Ide, A. Kashihara, M. Kano, H. Mariko, R. Masui, N. Nak-agawa, S. Yokoyama, S. Kuramitsu, and A. Konagaya, “Inference of S-system models of genetic networks using a cooperative coevolutionary algorithm,” Bioinform., vol. 21, pp. 1154–1163, 2005.

[13] C. G. Moles, P. Mendes, and J. R. Banga, “Parameter estimation in biochemical pathways: A comparison of global optimization methods,” Genome Res., vol. 13, pp. 2467–2474, 2003.

[14] N. Noman and H. Iba, “Inference of gene regulatory networks using S-system and differential evolution,” in Proc. Conf. Genetic Evol. Comput., 2005, vol. 1, pp. 439–446.

[15] K. Y. Tsai and F. S. Wang, “Evolutionary optimization with data colloca-tion for reverse engineering of biological networks,” Bioinform., vol. 21, pp. 1180–1188, 2005.

[16] S. J. Wu, C. T. Wu, and T. T. Lee, “Computation intelligent for eukaryotic cell-cycle gene network,” in Proc. IEEE Eng. Med. Biol. Soc. Conf., 2006, pp. 2017–2020.

[17] P. K. Liu and F. S. Wang, “Hybrid differential evolution with geometric mean mutation in parameter estimation of bioreaction systems with large parameter search space,” Comput. Chem. Eng., vol. 33, pp. 1851–1860, 2009.

[18] F. S. Wang and P. K. Liu, “Inverse problems of biochemical systems using hybrid differential evolution and data collocation,” Int. J. Syst. Synthetic Biol., vol. 1, pp. 21–38, 2010.

[19] S. Kikuchi, D. Tominaga, M. Arita, K. Takahashi, and M. Tomita, “Dy-namic modeling of genetic networks using genetic algorithm and S-system,” Bioinform., vol. 19, pp. 643–650, 2003.

[20] E. O. Voit and J. Almeida, “Decoupling dynamical systems for pathway identification from metabolic profiles,” Bioinformatics, vol. 20, pp. 1670– 1681, 2004.

[21] S. Y. Ho, C. H. Hsieh, F. C. Yu, and H. L. Huang, “An intelligent two-stage evolutionary algorithm for dynamic pathway identification from gene expression profiles,” IEEE/ACM Trans. Comput. Biol. Bioinform., vol. 4, no. 4, pp. 648–660, Oct./Dec. 2007.

[22] O. R. Gonzalez, C. K¨uper, K. Jung, P. C. Naval, E. Mendoza, Jr., and E. Mendoza, “Parameter estimation using simulated annealing for

S-system models of biochemical networks,” Bioinform., vol. 23, pp. 480– 486, 2007.

[23] C. M. Chen, C. Lee, C. L. Chuang, C. C. Wang, and G. S. Shieh, “Inferring genetic interactions via a nonlinear model and an optimization algorithm,” BMC Syst. Biol., vol. 4, no. 16, 2010.

[24] Y. Matsubara, S. Kikuchi, M. Sugimoto, and M. Tomita, “Parameter es-timation for stiff equations of biosystems using radial basis function net-works,” BMC Bioinform., vol. 7, no. 230, 2006.

[25] H. Murata, M. Koshino, M. Mitamura, and H. Kimura, “Inference of S-system models of genetic networks using product unit neural networks,” in Proc. IEEE Conf. Syst. Man Cybern., 2008, pp. 1390–1395.

[26] R. Xu, D. C. Wunsch II, and R. L. Frank, “Inference of genetic regulatory networks with recurrent neural network models using particle swarm opti-mization,” IEEE Trans. Comput. Biol. Bioinform., vol. 4, no. 4, pp. 1545– 5963, Oct./Dec. 2007.

[27] D. Thieffry, A. M. Huerta, E. Perez-Rueda, and J. Collado-Vides, “From specific gene regulation to genomic networks: A global analysis of tran-scriptional regulation in escherichia coli,” BioEssays, vol. 20, pp. 433– 440, 1998.

[28] P. K. Liu and F. S. Wang, “Inference of biochemical network models in S-system using multiobjective optimization approach,” Bioinform., vol. 24, pp. 1085–1092, 2008.

[29] D. S. Falconer and T. F. C. MacKay, Introduction to Quantitative Genetics, 4th ed. White Plains, NY: Longman, 1996.

[30] C. L. Ko, E. O. Voit, and F. S. Wang, “Estimating parameters for general-ized mass action models with connectivity information,” BMC Bioinform., vol. 10, no. 140, 2009.

[31] I. C. Chou and E. O. Voit, “Recent developments in parameter estimation and structure identification of biochemical and genomic systems,” Math. Biosci., vol. 219, pp. 57–83, 2009.

[32] W. S. Hlavacek and M. A. Savageau, “Rules for coupled expression of regulator and effector genes in inducible circuits,” J. Mol. Biol., vol. 255, pp. 121–139, 1996.

[33] M. J. Gacto, R. Alcala, and F. Herrera, “Integration of an index to preserve the semantic interpretability in the multiobjective evolutionary rule selec-tion and tuning of linguistic fuzzy systems,” IEEE Trans. Fuzzy Syst., vol. 18, no. 3, pp. 515–531, Jun. 2010.

[34] C. J. Carmona, P. Gonzalez, M. J. del Jesus, and F. Herrera, “NMEEF-SD: Non-dominated multiobjective evolutionary algorithm for extracting fuzzy rules in subgroup discovery,” IEEE Trans. Fuzzy Syst., vol. 18, no. 5, pp. 958–970, Oct. 2010.

[35] G. P. Rangaiah, Multi-Objective Optimization: Techniques and Applica-tions in Chemical Engineering. Hackensack, NJ: World Scientific, 2009. [36] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: NSGA-II,” IEEE Trans. Evol. Comput., vol. 6, no. 2, pp. 182–197, Apr. 2002.

[37] K. Deb, “Current trends in evolutionary multi-objective optimization,” Int. J. Simul. Multidisci. Design Optim., vol. 1, pp. 1–8, 2007.

[38] M. Sakawa and K. Kato, “An interactive fuzzy satisficing method for mul-tiobjective nonlinear integer programming problems with block-angular structures through genetic algorithms with decomposition procedures,” Adv. Oper. Res., vol. 2009, 2009.

[39] K. Deb, A. Sinha, P. J. Korhonen, and J. Wallenius, “An interactive evo-lutionary multiobjective optimization method based on progressively ap-proximated value functions,” IEEE Trans. Evol. Comput., vol. 14, no. 5, pp. 723–739, Oct. 2010.

[40] A. Sukstrienwong, “Solving multiobjective optimization under bounds by genetic algorithms,” Int. J. Comput., vol. 5, pp. 18–25, 2011.

[41] P. Pulkkinen and H. Koivisto, “A dynamically constrained multiobjective genetic fuzzy system for regression problems,” IEEE Trans. Fuzzy Syst., vol. 18, no. 1, pp. 161–177, Feb. 2010.

[42] S. J. Wu, C. H. Chou, C. T. Wu, and T. T. Lee, “Inference of genetic network of xenopus frog egg: Improved genetic algorithm,” in Proc. IEEE Eng. Med. Biol. Soc. Conf., 2006, pp. 4147–4150.

Shinq-Jen Wu received the B.S. degree in chem-ical engineering from National Taiwan University, Taipei, Taiwan, in 1986, the M.S. degree in chemi-cal engineering from National Tsing-Hua University, Hsinchu, Taiwan, in 1989, the M.S. degree in electri-cal engineering from the University of California, Los Angeles, in 1994, and the Ph.D. degree in electrical and control engineering from National Chiao-Tung University, Hsinchu, in 2000.

From September 1989 to July 1990, she was with the Laboratory for Simulation and Control

intelligent modeling techniques.

Dr. Wu is a member of the Phi Tau Phi Scholastic Honor Society. She is the Editor of Advances in Fuzzy Sets and Systems (Pushpa, Allahabad, India). Her name is included in Asian Admirable Achievers, Asia/American Who’s Who, Asia/Pacifica Who’s Who, and in Marquis Who’s Who in Science and Engineer-ing/in the World/in America/in Asia. She is a Life Fellow of the International Biographical Association. She is a Scientific Adviser to the International Bio-graphical Centre (IBC) Director General. She received the 21st Century Award for Achievement from the IBC and The Albert Einstein Award of Excellent from the American Biographical Institute, Inc.

Cheng-Tao Wu received the M.S. degree in electri-cal engineering from Da-Yeh University, Changhua, Taiwan, in 2006.

He is currently with the Department of Electrical and Control Engineering National Chiao-Tung Uni-versity, Hsinchu, Taiwan. He is involved in projects concerning optimization and detection in control systems.

Electrical and Control Engineering, NCTU, where he is currently a Professor. His current research interests include neural fuzzy systems, video processing and surveillance, and bioinformatics.