
\documentclass[conference]{IEEEtran}
% If the IEEEtran.cls has not been installed into the LaTeX system files, 
% manually specify the path to it: e.g.,
% \documentclass[conference]{../sty/IEEEtran} 

\usepackage{graphicx,times,psfig,amsmath} % Add all your packages here
\graphicspath{{./images/}}


% correct bad hyphenation here
\hyphenation{op-tical net-works semi-conduc-tor IEEEtran}

\IEEEoverridecommandlockouts	% to create the author's affliation portion
				% using \thanks
				
\textwidth 178mm   	% <------ These are the adjustments we made 10/18/2005
\textheight 239mm	% You may or may not need to adjust these numbes again
\oddsidemargin -7mm
\evensidemargin -7mm
\topmargin -6mm 
\columnsep 5mm

\begin{document}


% paper title: Must keep \ \\ \LARGE\bf in it to leave enough margin.
\title{\ \\ \LARGE\bf Bio-inspired Reverse Engineering of Regulatory Networks\thanks{CRAB Lab,
Department of Computer and Information Science (IDI),
Norwegian University of Science and Technology (NTNU), Norway; emails: csantini@idi.ntnu.no, pauline@idi.ntnu.no and gunnar.tufte@idi.ntnu.no).}
}

\author{Cristina Costa Santini, Gunnar Tufte and Pauline Haddow}
% avoiding spaces at the end of the author lines is not a problem with
% conference papers because we don't use \thanks or \IEEEmembership
% use only for invited papers
%\specialpapernotice{(Invited Paper)}

% make the title area
\maketitle

\begin{abstract}
Regulatory networks are complex networks.
In this paper the challenge of modelling these complex networks is addressed. 
A modelling architecture based on bio-inspired techniques applied to the search for Boolean network representations of sought regulatory networks is presented.
Two bio-inspired techniques are investigated when applied to the reverse engineering of a Boolean network model: a Genetic Algorithm and an indirect reverse engineering method. 
The latter is proposed as a means of addressing the challenge of modelling large and complex networks.
\end{abstract}

% no key words

\section{Introduction}\label{intro}

Modelling regulatory networks is an ongoing challenge in the field of systems biology due to the inherent complexity of such networks. Regulatory networks consist of both
signal transduction networks and genetic transcription networks. Signal transduction networks transmit signals from the boundary of a cell to the nucleus in a combinatorial process and genetic transcription networks control if and how genes are transcribed.

The complexity inherent in signal transduction pathways arises, not only from the size i.e. the number of pathways and components, but from the crosstalk between pathways and from the convergence of pathways. Further, particular genes influence multiple functions (pleitropy), providing a further level of complexity within regulatory networks~\cite{huang04}. Such complexity was confirmed by the publication of the human genome in 2001. Much fewer genes than originally estimated were shown to be present thus illustrating the fact that genes must  interact in even more complex ways in order to maintain homeostasis~\cite{venter01}.

Any form of modelling of biological processes needs to address the abstraction level sought i.e. the level of detail to be included and the scope or aim of the model~\cite{kitano02}. Further, modelling the dynamics inherent in complex biological processes requires vast resources. Although unlimited resources may be assumed, realistically effective computation platforms are required.

A further challenge inherent in modelling biological processes is the information available to modellers. On one side, there may be a vast amount of information or unstructured data, providing the modeller with a challenge to extract the relevant data. On the other side, data may not be complete thus providing a challenge as to how to cope with incomplete information.

In general, modelling approaches to regulatory networks may be said to attempt to provide a one to one description of the network or sub-network under investigation. As the complexity and the size of the network increases, challenges with respect to both the model creation and the resources required to run the dynamics of the model increase.

In this paper, an indirect reverse engineering method is proposed to address the challenges highlighted.
It is investigated through the reverse-engineering of a simple model, the yeast cell-cycle network model~\cite{li04bio}.
%This work extends and refines Haddow's approach~\cite{haddow08} where
%bio-inspired techniques are proposed to search for large-scale regulatory network
%models and the dynamics of such models are proposed simulated on a highly parallel hardware platform. 


This paper is organised as follows:
Section~\ref{motivation} presents the motivation and the context of the work.
Section~\ref{overview} presents an overview of the modelling architecture. 
Section~\ref{techniques} presents the search techniques applied and Section~\ref{brn} presents the representation chosen for the regulatory networks: Boolean networks. 
In Section~\ref{results} experiments and results of the reverse engineering of the yeast cell-cycle network model are presented and discussed. 
Section~\ref{future} discusses the future work and Section~\ref{conclusion} concludes the paper.



\section{Motivation}
\label{motivation}

Reverse engineering is an important first step in the modelling of biological regulatory networks, when, even though microarray data is available, there is no mechanistic understanding of the regulatory processes that are triggered by a certain stimulus.


Approaches to modelling regulatory networks have been classified by Huang~\cite{huang04} as local (detailed): aiming to achieve a detailed description of the processes modelled and global (abstract): aiming at understanding the general principles of the processes modelled. The former approach enables computational experiments to be executed, providing a faster and cheaper experimental platform than biological experiments. On the other hand, a systems understanding of regulatory networks will only be achieved if the models are more abstract i.e. the latter approach, so that biological principles can be extracted~\cite{huang01,maini01}. A further aspect supporting the abstract approach is, as stated in Section~\ref{intro}, the lack of available and suitable data. Quantitative measurements of biological components in space and time are still insufficient to provide a fully determined network.

There are several techniques being applied to modelling biological networks following both approaches, as reviewed in~\cite{fisher07}. The more detailed, descriptive models can be seen as important `bridges' between the hypotheses generated by the higher level abstract modelling techniques and the actual experimental validation. The approach required depends on the purpose of the model building and on the scope and abstraction level of the problem~\cite{kitano02}. In this work, the abstract approach is chosen.

\section{Overview of the modelling architecture}\label{overview}

The goal is to achieve an understanding  of the biological principles that govern regulatory networks. 
Figure~\ref{fig:devarchitecture} shows the modelling architecture.
This extends the work of Haddow~\cite{haddow08} who proposed a bio-inspired approach for searching for suitable regulatory network structures where the dynamics of such sought networks are studied on a highly parallel hardware platform. The goal being to generate realistic hypotheses by a circular process involving several model attempts, hypotheses generation and biological feedback. 
\begin{figure} [!tbp]
\begin{center}
\includegraphics[width=0.48\textwidth]{devarchitecture}
\caption{Overview of the modelling architecture.}
\label{fig:devarchitecture}
\end{center}
\end{figure}

There are a number of unique elements in the approach proposed by Haddow~\cite{haddow08}. 
It suggests evolutionary algorithm techniques (see section~\ref{EA}) to search for network solutions and introduces a newer technique termed artificial development (see section~\ref{artDev}) so as to scale up the search to large-scale network solutions. A further unique element is the computational platform itself. To meet the processing challenges of simulating the dynamics of large-scale networks, a highly parallel hardware platform is proposed. 

The modelling architecture is refined herein and, in particular, Boolean networks are introduced as the representation of the networks sought and an indirect reverse engineering method is proposed and investigated through the reverse-engineering of a simple model, the yeast cell-cycle network model~\cite{li04bio}.
This is a small network, allowing this preliminary reverse engineering experiment to be implemented in software.

Nevertheless, the final aim of the work is to reverse engineer large-scale networks that are going to be executed in hardware.
Therefore, a criteria the representation to be chosen had to fulfill was to be suitable for hardware implementation. Being based on simple processing elements, Boolean networks are suitable for a hardware implementation (see Section~\ref{future}).

In this work, the two different bio-inspired techniques, evolutionary algorithms and artificial development, are applied as a search engine to find possible Boolean regulatory networks (BRN) (see Section~\ref{brn}).
These are executed and  the state space analysis provides, for the case study in hand, an evaluation of their dynamics, i.e. a `fitness' measure. Such a process is a cyclic process, as shown in Figure~\ref{fig:devarchitecture}, stopping when some search criteria is met.  The output being a BRN hypothesis. Refinement of the hypotheses is again a cyclic process where biologists are involved to evaluate the hypotheses and provide revised input to the modelling process, as shown. The following sections provide more details about the elements of the modelling architecture.

\section{Bio-inspired techniques}\label{techniques}

\subsection{Evolutionary Algorithms (EA)}\label{EA}

Taking inspiration from Darwin's theories of evolution, evolutionary algorithms (EAs)~\cite{HOL75} aim to exploit the power of biological evolution as a computational search and/or optimisation tool. The artificial evolution process (evolutionary algorithms) is a dynamic process where at a given point in time one has what is termed in biology a generation: a population of individuals for the given species. Biologically, a generation change is not an event but a constantly unfolding process. However, in artificial evolution, a new generation is introduced through the application of reproduction and mutation techniques to selected individuals of the current generation. Further, unlike biological evolution there is a specified goal i.e.\ achieving a correct `fit' solution to the problem in hand.

EAs are challenged to produce scalable solutions due to their resource-greedy nature~\cite{DROSTEACI03}.  One source of the scalability problems of EAs is the use of a one-to-one mapping from genotype to phenotype. That is, the solution generated by the EAs may be directly mapped to the candidate solution. This means that all properties of the phenotype will need to be directly encoded in the genotype. A large complex phenotype will, therefore, require a large complex genotype.

\subsection{Artificial Development (AD)}\label{artDev}


Nature's way of handling complexity clearly points in the direction of a non one-to-one mapping from genotype to phenotype. Biological development is nature's way of coping with scaling. It provides a plan of the phenotype, describing how the organism is to be built rather than a blueprint containing explicit information about every detail of the
proposed phenotype. 


Unlike Evolutionary Algorithms there are no established methods for Artificial Development that may be applied in this work. 
Many efforts have been made in the field to test out different approaches to AD. However, little consensus has been achieved as to how such an approach should be applied and, in particular, little knowledge has been gained as to how different application areas should be approached. As such, a key challenge addressed herein is to find a way to apply Artificial Development to the reverse engineering of the Boolean network representation of regulatory networks.

Therefore, in order to address this challenge, a method is proposed and investigated, and is referred to as `Artificial Development reverse engineering method'. 
In this method, the genotype is updated (or tuned) according to the phenotype's performance. More details about this method are given in Section~\ref{resultsb}, where we use the reverse engineering $in-silico$ experiments as an example.


\section{Boolean Regulatory Networks}\label{brn}

Boolean networks can be used as a representation for regulatory networks.
A thorough discussion can be found on the seminal papers on the subject as well as in recent reviews and perspectives
~\cite{kauffman69,kauffman93,huang99,huang00,huang00b,huang01,huang04,shmulevich02,bornholdt05}.


\begin{figure} [!tbp]
\begin{center}
\includegraphics[width=0.48\textwidth]{architecture}
\caption{(a) Basic structure of a Boolean regulatory network.
(b) Scheme of the state space with the three attractors, illustrating the `attractor landscape' (taken from~\cite{huang99}) in the notation proposed by~\cite{wuensche98}.
}
\label{fig:architecture}
\end{center}
\end{figure}

Boolean networks were first introduced by Kauffman~\cite{kauffman69}.
In the Boolean idealization a regulatory network consists of a number \textit{n} of elements representing genes/proteins, which can be either active `on' (1) or inactive `off' (0).
The basic structure of a Boolean regulatory network is illustrated in Figure~\ref{fig:architecture}a.
Each square represents an element of the network.
These elements are interconnected as defined by the wiring diagram.
Each element receives \textit{k} inputs.
A Boolean function is assigned to each element, dictating how the input is processed into an output.
In the example shown the element C has two inputs ($k$=2) and its state is updated according to the $OR$ Boolean function.
The network state $S(t)$ is defined by the values of all the \textit{n} elements of the network at a given time \textit{t}.
>From the example, note that C=0 at time t and becomes 1 at time t+1 due to the values of A and B at time t. Their new values at t+1 reflect their own updates. 
By executing the network, i.e. executing the Boolean rules, a state transition occurs and all the states in the network are updated based on such transition rules such as the OR rule in the example. 

%The dynamics of this Boolean network will determine the various cell types~\cite{kauffman93} or cell states/fates~\cite{huang99} (attractor states of the network) such as differentiation, proliferation and apoptosis in an organism. 
%(Huang's view of cell \emph{fates} as attractors is complementary to Kauffman's proposal that cell \emph{types} are attractors, assuming that a cell type comprises a set of `neighbouring' attractors in which each attractor represents a cellular state~\cite{huang00b}.)
The dynamics of this Boolean regulatory network will determine the various cell fates (attractor states of the network) such as differentiation, proliferation and apoptosis~\cite{huang99}.
Figure~\ref{fig:architecture}b shows the scheme of the state space with the three attractors.
Every grey dot in the state space represents a unique, transient (i.e., unstable) network state at a certain time (illustrated by the two external arrows on the left). These network states are connected
by arrows that form the trajectories. The state
space is divided (dotted line) into three basins of attraction with an
attractor (black dots) at the centre of each. The size of the basins of attractions is correlated with the
stability of the respective attractors. The arrows at the attractor
boundaries and their size illustrate the probabilities of transition
between the attractors.

%\subsection{Reverse engineering Boolean regulatory networks}

Examples can be found in the literature that experimentally support this hypothesis of cell fates as attractor states. In addition, they provide further evidence that the Boolean network representation, even though abstract and simplistic, captures the essential aspects of regulatory networks.
In~\cite{huang00}, Huang and Ingber simulate the dynamics of the signalling system within capillary endothelial cells as a dynamic Boolean network and show that the attractors correspond to the cell's fates.
In the work shown in~\cite{li04bio}, the authors propose and investigate a dynamic Boolean network model for the yeast cell-cycle (see Section~\ref{results} for more details on this model).
A further example of this approach is presented in ~\cite{helikar07}, where the authors create a discrete Boolean model of signal transduction and they test the model's ability to replicate known qualitative behaviour. They support the idea that signal transduction networks are non-trivial decision making systems capable of clustering widely varying input combinations into a few cell states.



\section{Experiments and Results}
\label{results}

In order to investigate the bio-inspired techniques when applied to the reverse engineering of regulatory networks, we performed preliminary $in-silico$ experiments aiming at the reverse engineering of the yeast cell-cycle network model, presented in~\cite{li04bio} and reproduced in Figure~\ref{fig:yeastnet1}. 
This model was chosen for being a simple BN model that, even though simple, models successfully the yeast cell-cycle phases and stationary states.

\begin{figure} [!tbp]
\begin{center}
\includegraphics[width=0.48\textwidth]{yeastnet1}
\caption{Yeast cell-cycle network. From~\cite{li04bio}.
}
\label{fig:yeastnet1}
\end{center}
\end{figure}

In the model, each node \textit{i} has only two states, \textit{S$_{i}$}=1 and \textsl{S$_{i}$}=0, representing the active and inactive state of the protein, respectively.
The protein states in the next time step are determined by the protein states in the present time step via the following transition rule:

\begin{equation}
  \label{eq:statesupdate}
  S_i(t+1) = \left\{ \begin{array}{ll}
      1,              & \textrm{$\sum_{j}a_{ij}S_{j}(t) > 0$}\\
      0,						  & \textrm{$\sum_{j}a_{ij}S_{j}(t) < 0$}\\
      S_{i}(t),						  & \textrm{$\sum_{j}a_{ij}S_{j}(t) = 0$}
    \end{array} \right.
\end{equation}



For the results presented here \textit{a$_{ij}$}=1 for a green arrow from protein \textit{j} to protein \textit{i} and \textit{a$_{ij}$}=-1 for a red arrow from \textit{j} to \textit{i}. 
A green arrow represents positive regulation and a red arrow represents `deactivation' (inhibition, repression, or degradation).

The model also has `self-degradation' (yellow loops) on those nodes that are not negatively regulated by others. This self-degradation is modeled as a time-delayed interaction: if a protein with a self yellow arrow is active at time t and if its total input is zero at t+1, it will be degraded, i.e. \textsl{S$_{i}$(t+1)}=0.
Using this dynamic model, the attractors of the network can be found; each of the 2$^{11}$=2048 initial states in the 11-node network eventually flow into one of seven stationary states (fixed points). Those are shown in Table 1, where each fixed point is represented in a row. The first column is the size of the basin of attraction for the fixed point; the other 11 columns show the protein states of the fixed point\footnote{Please refer to~\cite{li04bio} for more information about the model and its dynamic properties, as well as for the results on network perturbations.
It can be mentioned that we reproduced and extended with success those experiments.
}.


\begin{table*}[!tbp]\footnotesize
\caption{The fixed points of the cell-cycle network. From \cite{li04bio}}
\begin{center}
\begin{tabular*}{0.95\textwidth}{@{\extracolsep{\fill}}@{}c@{}@{}c@{}@{}c@{}@{}c@{}@{}c@{}@{}c@{}@{}c@{}@{}c@{}@{}c@{}@{}c@{}@{}c@{}@{}c@{}}
\multicolumn{1}{c}{Basin size}
& \multicolumn{1}{c}{Cln3}% time}
& \multicolumn{1}{c}{MBF}% time}
& \multicolumn{1}{c}{SBF}% cost}
& \multicolumn{1}{c}{Cln1,2}% cost}
& \multicolumn{1}{c}{Cdh1} 
& \multicolumn{1}{c}{Swi5} 
& \multicolumn{1}{c}{Cdc20} 
& \multicolumn{1}{c}{Clb5,6} 
& \multicolumn{1}{c}{Sic1} 
& \multicolumn{1}{c}{Clb1,2} 
& \multicolumn{1}{c}{Mcm1} \\% cost}\\
\hline
\footnotesize
%& time&time& cost&cost&cost& time&time& cost&cost&cost&cost&cost\\ \hline
 $1,764$   &  $0$  &  $0$  &  $0$ & $0$ &  $1$  &  $0$  &  $0$ & $0$ &  $1$  &  $0$  &  $0$  \\ 
 $151$   &  $0$  &  $0$  &  $1$ & $1$ &  $0$  &  $0$  &  $0$ & $0$ &  $0$  &  $0$  &  $0$  \\ 
 $109$   &  $0$  &  $1$  &  $0$ & $0$ &  $1$  &  $0$  &  $0$ & $0$ &  $1$  &  $0$  &  $0$  \\ 
 $9$   &  $0$  &  $0$  &  $0$ & $0$ &  $0$  &  $0$  &  $0$ & $0$ &  $1$  &  $0$  &  $0$  \\ 
 $7$   &  $0$  &  $1$  &  $0$ & $0$ &  $0$  &  $0$  &  $0$ & $0$ &  $1$  &  $0$  &  $0$  \\ 
 $7$   &  $0$  &  $0$  &  $0$ & $0$ &  $0$  &  $0$  &  $0$ & $0$ &  $0$  &  $0$  &  $0$  \\ 
 $1$   &  $0$  &  $0$  &  $0$ & $0$ &  $1$  &  $0$  &  $0$ & $0$ &  $0$  &  $0$  &  $0$  \\ 
 \hline
\end{tabular*}
\label{table:basins}
\end{center}
\end{table*}



The aim of the experiments presented is to investigate the feasibility of the proposed AD reverse engineering method. 
Therefore, reverse engineering experiments when applying this method are shown.
In addition, in order to get an insight into the problem at hand, i.e.\ to investigate how hard the problem is, another bio-inspired method is also applied: a canonical Genetic Algorithm with standard parameters.


The data set for the $in-silico$ experiments include only partial data about the network under investigation. 
The choice of only including partial data is motivated by the fact that in reality, i.e.\ when investigating a biological system, it is unlikely that it will be possible to have inform about all the initial conditions and corresponding attractors. The full data set for the 11-node network consists of 2048 initial states. In the $in-silico$ experiments shown here, information about only 100 of the possible initial states is provided.

In addition, 
%even more relevant to a biological application, 
we show the results of an $in-silico$ experiment when applying the AD reverse engineering method and when the data set includes only a few trajectories (time series data). A relevant investigation to biological applications.

Because the GA is applied directly to the Boolean network representation of the yeast cell-cycle network presented in this Section, we first show the GA results, followed by the results for the AD reverse engineering method.








\subsection{Applying EAs}
\label{resultsa}

The Evolutionary Algorithm we used is a canonical Genetic Algorithm (GA); it has a one-to-one genotype-phenotype mapping and its genotype is evolved.

Figure~\ref{fig:yeastBNGA} shows the table for the Boolean network representation that corresponds to the wiring diagram of the yeast cell-cycle network shown in Figure~\ref{fig:yeastnet1}.
Every column represents a network element, or node, and lines 1 to 5 inform which are the inputs to that node, as well as if they are positive regulations or a `deactivations'.
This representation shown in Figure~\ref{fig:yeastBNGA} is what is called here in this experiment \textit{the genotype}.

%Lines 6-10 inform if it is a positive regulation (green arrows in Figure~\ref{fig:yeastnet1}) or if it represents `deactivation' (inhibition, repression, or degradation) (red arrows in Figure~\ref{fig:yeastnet1}). Line 11 inform if that node suffers from `self-degradation' or not (yellow loops in Figure~\ref{fig:yeastnet1}).
\begin{figure} [!tbp]
\begin{center}
\includegraphics[width=0.5\textwidth]{yeastBNGA}
\caption{Table showing the yeast cell-cycle network elements and their respective positive regulations or `deactivations'.}
\label{fig:yeastBNGA}
\end{center}
\end{figure}

In order to perform a reverse engineering investigation, we have to assume that some biological information about the regulatory network is known a priori, but not all. A possible way of doing that is by specifying an initial representation, i.e. an initial genotype.
For the results shown in this section, 
%in the Boolean network representation of the yeast cell-cycle network model, 
the inputs to nodes 2 to 7 are searched for, whereas the inputs to the other nodes are considered known biological information.

The GA will then search for the inputs to the unknown nodes.
Consider node 2 for example: during the GA search it can be determined that it is positively regulated by nodes 4 and 6 and negatively regulated by node 11.
Then the fitness of this network is determined by executing it to all and each of the given initial conditions and measuring the hamming distance between the obtained final attractor and the yeast network's final attractor for that initial condition, which can only be one of the seven attractors shown in Table~\ref{table:basins}.

The executed network, or the dynamic network, or in nonlinear dynamics terms, the state space, is what is called here in this experiment \textit{the phenotype}.

The average hamming distance for all the given initial conditions is then the phenotype distance \textit{d$_P$}, with the fitness being $F=e^{-d_P}$~\cite{fernandez07}.\\

\subsubsection{EAs Results}

The parameters of the \textit{in-silico} experiment\footnote{This \textit{in-silico} experiment is implemented in Matlab 7.6.0 and its execution time is in the hours time span.} are: 
\begin{description}
	\item[- generations: 1,000;]
	\item[- initial population: 2,000;] 
	\item[- elitism: 1,000;]
	\item[- mutation: 600;]
\end{description}
The first generation has 2,000 randomly created network solutions.
The following generations are formed by the 1,000 best networks of the previous generation, 600 mutated networks and 400 new randomly created networks.
Mutation is performed on 600 randomly chosen elements of the previous generation, not on the 600 best ones.
With those parameters, a maximum of 10$^6$ different proposed networks can be executed and evaluated in a run.



\begin{figure} [!tbp]
\begin{center}
\includegraphics[width=0.45\textwidth]{graphGAred}
\caption{Applying GA and 100 initial conditions. 
The green series shows the maximum fitness given 100 initial conditions and the blue series shows the fitness when the network found is evaluated for the whole state space, i.e. with 2048 initial conditions; the red series shows the number of networks evaluated, in each of the 10 runs.}
\label{fig:GAb}
\end{center}
\end{figure}

The results to 10 runs are graphically shown in Figure~\ref{fig:GAb}.
The green series shows the maximum fitness given 100 initial conditions, with average fitness 0.90. The exact solution for this data (fitness 1) was found in 30\% of the runs (runs 2, 5 and 9). The fitness of the networks found in each run when executed and evaluated for the whole state space, i.e. for the 2048 initial states, is shown in the blue series, and is on average 0.75. Also, on average, a solution was found around generation 800, i.e. after 800x1,000 different proposed networks have been executed.







\subsection{Applying the Artificial Development (AD) reverse engineering method}
\label{resultsb}


We propose and investigate an AD reverse engineering method in which the biological network is represented by a list of influences, interactions or regulations. 
Here in this experiment this list is called \textit{the genotype}. It carries the information about all the possible influences in the system,
%Following this approach, the genotype carries the information about all the possible influences in the system, 
and informs the probability in which that influence is going to be present in the list of influences that uniquely represent the network, which for now on will be referred to as \textit{the network}.
For example, consider Figure~\ref{fig:yeastnet1}, in which each arrow represents an influence, an activation or a deactivation: \textit{`MBF (node 2) is positively regulated (green arrow) by Cln3 (node1)'} is one of the influences that compose the network.

\begin{figure*} [!tbp]\footnotesize
%\begin{center}
(a) Initial genotype with the probabilities for each influence $i$\\
\includegraphics[width=0.98\textwidth]{genophenotypea}\\
(b) A possible created network with the number of occurrences and with the normalized presence of each influence $i$\\
\includegraphics[width=1\textwidth]{genophenotypeb}\\
(c) Final genotype for run 5 with the updated probabilities for each influence $i$\\
\includegraphics[width=0.98\textwidth]{genophenotypec}
\caption{Artificial Development reverse engineering method. (a) Table explicitly showing the initial genotype: nodes in row 1 are being (positively or negatively) influenced by the nodes in row 2; row 3 informs the probability with which each influence will be present in the network; (b) Table showing a possible created network: row 3 informs the number of times influence $i$ appears in the network and row 4 is the normalization of row 3; (c) Table showing the final genotype for one of the runs: row 3 informs the final updated probabilities for each influence $i$.}
\label{fig:genophenotype}
%\end{center}
\end{figure*}

Figure~\ref{fig:genophenotype}a explicitly shows the genotype used in the experiments shown here. 
The influence mentioned above, for example, is the first influence in the genotype.

As the reverse engineering investigation has to mimic a real biological situation, we suppose we only have some information about the network. In the experiment shown here, this is done by assigning different propabilities to the influences in the genome and by including some random influences that are not in the network.
It can be observed that the first 29 influences in the genotype correspond to the interactions present in the yeast cell-cycle network (Figure~\ref{fig:yeastnet1}). The other influences (columns 30 to 39) are included as possible influences.
Consider the third row of the table in Figure~\ref{fig:genophenotype}a: the influences with probability 1 correspond to known information and will certainly be in the network and this probability will never change. The influences with probability less than 1 will have their probabilities updated according to the `performance' of the network.



The adaptation process is depicted in Figure~\ref{fig:adaptation} and proceeds as follows: given the genotype, a network is built; given the genotype shown in Figure~\ref{fig:genophenotype}a, a possible network is created. This is shown in Figure~\ref{fig:genophenotype}b. This process can follow diverse strategies. For the results shown here, all the genotype is checked twice, i.e. all the influences in the genotype have their respective probabilistic chances to be added to the network twice.
The outcome is shown in the third row of the table in Figure~\ref{fig:genophenotype}b which informs how many times the respective influence is present in the network.

\begin{figure} [!tbp]
\begin{center}
\includegraphics[width=0.45\textwidth]{adaptation}
\caption{The iterative process of the Artificial Development reverse engineering method proposed.}
\label{fig:adaptation}
\end{center}
\end{figure}

Then, just as with the experiments using GA, this network is executed according to (\ref{eq:statesupdate}) for each given initial condition. 
To execute the network means to investigate its dynamics for all the provided initial conditions. The result of this execution is the phenotype (in nonlinear dynamics terms, the state space). 
Following the method proposed here, \textit{the phenotype} is therefore the product between \textit{the genotype} (all possible influences and corresponding probabilities of being present in \textit{the network}) and \textit{the environment} (represented by the initial conditions).

The phenotype is evaluated by measuring its distance to the yeast cell-cycle network's phenotype. 
If this phenotype distance is smaller than all the phenotype distances previously measured, the genotype is updated.
The update is a weighted update and it takes into account the current probability in the genotype ($prob\_G$ - Figure~\ref{fig:genophenotype}a third row) and the normalized list showing how many times a certain influence $i$ appeared on the phenotype ($norm\_presence\_phenotype$ - Figure~\ref{fig:genophenotype}b forth row).
For the results shown here the weight of the phenotype is 0.2. Therefore, the updated probability in the genotype for a certain influence \textit{i} is $updated\_prob\_G_i = prob\_G_i*0.8 + norm\_presence\_phenotype_i*0.2$.

This process, illustrated in Figure~\ref{fig:adaptation}, of creation and execution of the network and evaluation of the phenotype, followed by the possible update of the genotype, 
is repeated for a number of iterations or until a certain phenotype distance is achieved. %(Figure~\ref{fig:genophenotype}c).

The idea behind this process is the following: a certain influence that is present in a `good' phenotype should have an increased probability of being present in the network; in the same way, an influence that is not present in a good phenotype should have its probability of being in the network decreased.

Figure~\ref{fig:genophenotype}c shows the final genotype for one of the runs. Note that the list of interactions is the same as in the initial genotype (Figure~\ref{fig:genophenotype}a), but with the probabilities, shown in the third line, updated. After many updates, as expected, most of the probabilities up to influence 29 were updated towards 1, whereas the probabilities for most of the influences that are not in the original network (columns 30 to 39) were updated towards 0.


\paragraph{The experimental setup}

In this experiment 100 randomly chosen initial conditions are given. 
Figure~\ref{fig:ADb} shows the results for 10 runs. The target phenotype, which is the corresponding yeast cell-cycle network's final attractor for each of those 100 initial conditions, was found in all the 10 runs (green series).
The fitness for the network found in each run when executed and evaluated for the whole state space, i.e. for the 2048 initial states, is shown in the blue series and is on average 0.94. Most of the runs evaluated less than 3,000 proposed networks before the target was achieved (red series in Figure~\ref{fig:ADb}).



\begin{figure} [!tbp]
\begin{center}
\includegraphics[width=0.45\textwidth]{graphADred}
\caption{Applying the AD reverse engineering method and 100 initial conditions.
The green series shows the maximum fitness given 100 initial conditions and the blue series shows the maximum fitness when the network found is evaluated for the whole state space, i.e. with 2048 initial conditions; the red series shows the number of networks evaluated, in each of the 10 runs.}
\label{fig:ADb}
\end{center}
\end{figure}

It is important to add that, when providing data about all the 2048 (and not only 100) initial conditions and corresponding final attractors, the network that corresponds exactly to the yeast cell-cycle network is found in all the 10 runs and in most of the runs less than 10,000 proposed networks were evaluated before the target was achieved.


\paragraph{Further Investigations}

\begin{figure} [!tbp]
\begin{center}
\includegraphics[width=0.45\textwidth]{statespace4traj}
\caption{State-space graph of the yeast cell-cycle network model, taken from~\cite{willadsen07}.
Each node in the graph represents a state of the system, and an edge between
nodes represent a dynamic transition between states, with a loop showing a fixed point. The seven fixed points of the network are shown in Table~\ref{table:basins}.
The arrows illustrate the 4 trajectories that were given to the reverse engineering algorithm.
}
\label{fig:statespace4traj}
\end{center}
\end{figure}

Motivated by the type of time series data\footnote{Microarray mRNA data from 15 minutes to 24 hours after gastrin stimulation~\cite{bruland08}.} that is used by our collaborating group, the Gastrin Systems Biology group, another set of $in-silico$ experiments was performed.
Instead of providing as data for the reverse engineering process the initial conditions and evaluating the corresponding final states, in these experiments the data is composed of only 4 trajectories, as graphically exemplified in Figure~\ref{fig:statespace4traj}: two trajectories from the biggest basin of attraction and one from each of the two second biggest. In this situation, all the states in the four trajectories of the proposed networks are evaluated.
In all the 10 runs, a network was found that would match exactly (fitness 1) the activation of all the nodes for the four given trajectories.
Generalizing one of these solutions for the whole state space, i.e. for all the 2048 states, its fitness is 0.86.
On average, 18,000 proposed networks were investigated in each run.

Another equivalent simulation was performed with 6 trajectories: the same four trajectories from before plus two from the biggest basin. Again, the networks proposed could match exactly all the trajectories, for all the 10 runs, and 14,000 proposed networks were investigated on average, in each run. In addition, when evaluating one of these solutions for the whole state space, i.e. for all the 2048 states, its fitness is 0.95. This result suggests that when more data is given, the method finds a solution that better represents the whole state space.

Both results show that the method is applicable to time series data.


\subsection{Discussion}

The search space of the Boolean network increases exponentially with the number of nodes: $(n^k)^n$.
This implies that a sparsely connected network ($k$=2), with $n$=100 nodes, gives a total of $10^{400}$ possible networks. 

In the experiments shown here the indirect method applied managed to generate valid networks by only evaluating a fraction of all the possible networks in the total search space. 
As such, an indirect method would be favorable if large-scale networks consisting of hundreds or thousands of nodes are to be the target.

%The search space of the Boolean network increases exponentially with the number of nodes: $(n^k)^n$.
%As an example, to find the appropriate wiring diagram of a Boolean network with $n$=100 and low connectivity $k$=2, requires searching through $10^{300}$ possible networks.
%Therefore an indirect technique would be more suitable to search through the possible networks, as for large-scale networks, i.e. when $n$ is in the scale of hundreds or thousands, the search space gets considerably large.
%
Summarizing the results shown here, the indirect method, i.e. the Artificial Development reverse engineering method, investigated 3,000 networks and the best network had a fitness of 0.94. The Genetic Algorithm investigated 800,000 networks and the best network had a fitness of 0.75. 
Even though a direct comparison can't be made as neither method was exhaustively investigated in terms of their parameters, it can be noted that the number of networks investigated differ in two orders of magnitude. 
Therefore, the Artificial Development reverse engineering method shows promise but still needs to be investigated when applied to large-scale networks.

\section{Future Work}
\label{future}

In order to investigate the AD reverse engineering method proposed here when applied to large-scale networks, an artificial large-scale network~\cite{mendes03,bulcke06} can be used. Artificial networks were previously used to investigate other reverse engineering methods.
%~\cite{soranzo07
%,huangzan07
We plan at implementing this large-scale artificial network in hardware, as the dynamic execution of the network (or of some networks) in parallel will certainly represent a considerable speed-up of the reverse engineering process.

If a model of a Boolean network is to run on a sequential machine i.e.\ a standard microprocessor, as in the case of the results shown here, the run-time is given by the sequential update of the given number of nodes ($n$) in the network modelled. On the other hand, if the Boolean network is executed on a massive parallel hardware architecture, all nodes execute their Boolean function in parallel (one state change) and thus the run-time is not influenced by the number of nodes in the network. The simple processing element required by a Boolean network is only a simple Boolean function. 

Reconfigurable technology is an existing technology providing a suitable architecture with large amounts of simple processing units. Further the technology is flexible, enabling both many reconfigurations, a requirement of the bio-inspired techniques, and dynamic reconfiguration, allowing for updates in the hardware during simulation of the dynamics of the modelled network.  


\section{Conclusion}\label{conclusion}

This paper addresses the challenge of reverse engineering and modelling regulatory networks.
The Boolean network representation is supported as a suitable abstract representation that even though simplistic, would correspond to the modelling purposes.
An indirect reverse engineering method is proposed as a means of addressing the challenge of modelling large and complex networks.

%An important aspect of this architecture is the simple, straightforward way in which the modelling process sits into an iterative process between computer scientists and biologists, enabling the generated hypotheses to contribute to further refinements of the model.

A preliminary reverse engineering investigation is made through a set of $in-silico$ experiments, using two bio-inspired techniques: GA and the AD reverse engineering method proposed. Results show that both techniques are successful in reverse engineering the yeast cell-cycle network. However, on average, a better solution was found and less proposed networks were investigated when applying the AD reverse engineering method than when applying GA.






% conference papers do not normally have an appendix

% use section* for appendix
%\section*{Appendix}
%Put your appendix here if you have any.

% use section* for acknowledgment
\section*{Acknowledgment}
% optional entry into table of contents (if used)
%\addcontentsline{toc}{section}{Acknowledgment}
This work was carried out during the tenure of an ERCIM ``Alain Bensoussan'' Fellowship Programme.
C.C.S. thanks Prof. Hidde de Jong's group for profitable discussions during the ERCIM Exchange Programme.
We thank the Gastrin Systems Biology group at the Department of Cancer Research and Molecular Medicine, NTNU, for providing the biological realism.
 
% trigger a \newpage just before the given reference
% number - used to balance the columns on the last page
% adjust value as needed - may need to be readjusted if
% the document is modified later
%\IEEEtriggeratref{8}
% The "triggered" command can be changed if desired:
%\IEEEtriggercmd{\enlargethispage{-5in}}

% references section
% NOTE: BibTeX documentation can be easily obtained at:
% http://www.ctan.org/tex-archive/biblio/bibtex/contrib/doc/

% can use a bibliography generated by BibTeX as a .bbl file
% standard IEEE bibliography style from:
% http://www.ctan.org/tex-archive/macros/latex/contrib/supported/IEEEtran/testflow/bibtex
\bibliographystyle{IEEEtran}
% argument is your BibTeX string definitions and bibliography database(s)
%\bibliography{IEEEabrv,../bib/paper}
\bibliography{c:/ccsreferences}


\end{document}
