\chapter{Modelling Chemical Reaction Networks}
\label{modellingCSNs}
In Section \ref{sec:benefits} we outlined the potential applications of realising and evolving CIPNs \textit{in silico} as information processing devices. In order to evolve such CIPNs it is necessary to specify, represent, simulate and analyse these chemical reaction networks through the use of modelling techniques. In this chapter a review of different existing techniques for modelling chemical reaction networks is conducted. The outcome of this review is to select a suitable approach to model CIPNs in our project. This multi-disciplinary review results from the collaborative effort that we initiated with the bio-analysis group at the Friedrich Schiller University of Jena (UJ) in Germany (an ESIGNET partner). The survey on Markov chains, chemical master equations and SBML/CellML was conducted by the UJ group (which includes Dr. Thomas Hinze, Thorsten Lenser and Dr. Peter Dittrich). Bayesian networks, term rewriting systems, Petri nets, cellulat and agent-based/learning classifier systems techniques were reviewed by myself. Differential equations and $\pi$-calculus were examined by both the UJ group and myself. The evaluation criteria and comparison table were equally devised and realised by the UJ and Dublin City University based group (including Dr. George G. Mitchell and myself). This chapter combines some of the materials published in several ESIGNET deliverables and at various international conferences \citep{esignet31, esignet41, decraene2006msm, DBLP:conf/eurogp/DecraeneMMK07, decraene2008mae}.

\section{Introduction}
A variety of modelling techniques for biological reaction networks have been established in recent years \citep{alon2007isb}. We identify several main branches of modelling techniques:

\begin{itemize}

	\item \textit{Deterministic}: Chemical reactions are approximated as continuous deterministic processes at the macroscopic/system level. The system's variable states are uniquely determined by the pre-specified parameters describing the reactions (e.g., molecular concentration, reaction rates, etc.) and initial states of these variables. Given an initial set of pre-specified parameters, deterministic models enable one to monitor, predict and describe the dynamics of the system over time and/or space. Examples of deterministic modelling techniques include: ordinary/partial differential equations  \citep{hand,1015.34001,pmid15207904,pmid8816754}, Michaelis-Menten models \citep{heinrich1996rcs} and power-law models \citep{vera2007plm}.
	
	\item \textit{Stochastic}: In contrast with deterministic approaches, stochastic models explicitly account for the uncertainty that is involved in molecular processes. The system's variable states are determined by the pre-specified system's parameters and through the use of random variables. By addressing randomness or variability, stochastic models provide a more detailed representation of the system's potential dynamics (and not only the average behaviour as in deterministic approaches). Multiple executions of a stochastic model generate unique (from one another) dynamics/observations. The latter can be used to estimate probability distributions of the system's potential states (assisting in the construction of probabilistic models, see below).	Examples of stochastic modelling techniques include: Markov chains \citep{gomez2001ppu} and chemical master equations \citep{gillespie-01}.
	
	\item \textit{Probabilistic}: Here, the description of stochastic processes/data is addressed in terms of probability. Probabilistic modelling techniques are deterministic approaches which may infer probabilistic relationships between molecular species/system's states from empirical observations. In contrast with stochastic approaches, a probabilistic model is a statistical inference and description technique which does not represent the underlying stochastic molecular mechanics. Given the initial states of the molecular species, these approaches provide a probability-based description of the system's states. The predictive power of these techniques relies on the probabilistic distributions inferred by the model upon a range of \textit{in vivo/silico} experimental observations (i.e., the training set).	An example of probability modelling technique include: Bayesian networks \citep{sachs02} and hidden Markov models \citep{goutsias2006hmm}.
	
	\item \textit{Algebraic}: Modelling discrete characteristics of chemical reaction networks is principally achieved with algebraic approaches. A common basic assumption for these approaches is a finite or recursive enumerable number of elementary objects. Each object is considered as the smallest unit that can be processed by the system model. In particular, a definition of objects determines the granularity and abstraction level of corresponding models (hierarchically composed of objects, classes of objects, and temporal interaction rules). Both biomolecules and processes can form these objects. Interaction between these objects is usually specified by a relationship between system configurations. The whole system description is based on discrete transitions. This allows structural and comparative analysis of both system composition and behaviour, independent of numerical simulation results. Examples of algebraic modelling techniques include: P-systems \citep{paun1998dcn,paun2006mst}, broadcast language \citep{holland1992ana}, Alchemy \citep{fontana1994aft}, Boolean networks \citep{genoud1999}, $\pi$-calculus \citep{DBLP:conf/psb/RegevSS01} and Petri nets \citep{reddy1993pnr}.
	
	\item \textit{Agent-based}: Agent-based models (ABMs) extend the algebraic framework by introducing richer features in the computational units (i.e., agents). ABMs are commonly implemented with Object-Oriented programming environments in which agents are instantiations of object classes. The latter is a collection of properties (e.g., size, location, concentration, etc.) and methods (e.g., move, die, react, etc.). Agent-based simulations typically involve a large number of molecular and/or cellular agents which are executed in a concurrent or pseudo-concurrent manner. Each agent possesses its own distinct state variables, can be dynamically created/deleted and is capable of interacting with the other agents. The agents' computational methods may include stochastic processes resulting in a stochastic behaviour at the system level. Examples of agent-based modelling techniques include: Stochsim \citep{lenovere-shimizu-01}, Cellulat \citep{gonzalez2003cab} and AgentCell \citep{emonet2005ads}. A review of agent-based techniques is given by \cite{chavali2008cep}.

\end{itemize}

Deterministic and stochastic approaches are the most frequently employed and studied approaches in the field, whereas the attention given to the use of probabilistic, algebraic and agent-based approaches is more recent but rapidly growing. 
 
Our review of modelling techniques is not exhaustive. We select and review a limited number of techniques which exemplify the above families of modelling approaches. Following on from this we propose a model comparison table. We finally relate this evaluation with the requirements of our project to select a suitable technique with regard to the modelling, simulating and analysing of organisationally closed CIPNs. 

\section{Survey of modelling approaches}
\label{sec:surveyModel}
We review a selection of modelling techniques used in the study of CIPNs: differential equations, Markov chains, chemical master equations, Bayesian networks, Term Rewriting Systems, Petri nets, $\pi$-calculus, Cellulat and Agent-based Learning Classifier Systems. We then present the Systems Biology Markup Language (SBML) and CellML which allow one to specify and disseminate biochemical network models using a standardised language.  These markup languages also permit the migration of reaction network models between differing modelling approaches.

\subsection{Differential equations}
Chemical reactions are approximated as continuous deterministic processes at the macroscopic level. Differential equations provide a global understanding of a system and are commonly employed to model chemical reaction networks \citep{hand,1015.34001,pmid15207904,pmid8816754}. Given an initial set of pre-specified properties describing the reactions (e.g., molecular concentration, reaction rates, etc.), this modelling approach enables one to monitor, predict and describe the dynamics of the system over time and/or space. 

Here, state variables represent the concentrations of molecular species occurring in a well-stirred reactor with no in/out-flows. The following equation governs the dynamics of each species $S$  whose rate of change in concentration $[S]$ depends on the production and consumption rates $v_p$ and $v_c$:
\begin{equation}
\frac{d [S](t)}{d t}=v_{p}([S](t))-v_{c}([S](t)).
\end{equation}

In mass-action kinetics, these rates result from the reactant concentrations, their stoichiometric factors $a_{\{i,j\}} \in \mathbb{N}$ (reactants), $b_{i,j} \in \mathbb{N}$ (products) and kinetic constants $k_{j} \in \mathbb{R}_{+}$ assigned to each reaction quantifying its speed. For a reaction system with a total number of $n$ species and $r$ reactions

\begin{eqnarray*}
a_{11}S_1+a_{12}S_2+\ldots+a_{1n}S_n & \stackrel{k_1}{\longrightarrow} & b_{11}S_1+b_{12}S_2+\ldots+b_{1n}S_n\\
a_{21}S_1+a_{22}S_2+\ldots+a_{2n}S_n & \stackrel{k_2}{\longrightarrow} & b_{21}S_1+b_{22}S_2+\ldots+b_{2n}S_n\\
& \vdots & \\
a_{r1}S_1+a_{r2}S_2+\ldots+a_{rn}S_n & \stackrel{k_r}{\longrightarrow} & b_{r1}S_1+b_{r2}S_2+\ldots+b_{rn}S_n
\end{eqnarray*}
the corresponding ordinary differential equations (ODEs) read:
\begin{displaymath}
\frac{\mathrm{d}\, [S_i]}{\mathrm{d}\, t} = \sum\limits_{j=1}^{r} \left( k_{j}\cdot (b_{ji} - a_{ji}) \prod\limits_{h=1}^{n} [S_{h}]^{a_{jr}} \right)
\end{displaymath}
In order to obtain a concrete trajectory, all initial concentrations $[S_{i}](0) \in \mathbb{R}_{+}$, $i = 1, \ldots, n$ have to be specified. Solving this ODE system together with given initial values allows us to describe the temporal behaviour of the reaction system \citep{dittrich-articial}. 

Reaction-diffusion models take into account the spatial location of molecules and allow species concentrations in different spatial locations to vary continuously. These models are specified with sets of Partial Differential Equations (PDEs) \citep{fritz}. Solutions to PDEs derived from reaction-diffusion models provide an approximation of the species concentrations as a function $[S](t,x)$ of both time $t$ and space $x$:

\begin{equation}
\frac{\partial [S](t,x)}{\partial t} = D \frac{\partial ^2 [S](t,x)}{\partial x^2}-v([S](t,x))\frac{\partial [S](t,x)}{\partial x}+v_{p}([S](t,x))-v_{c}([S](t,x))
\label{eqPDE}
\end{equation}

Equation \ref{eqPDE} is an example PDE where the variables and functions represent: $[S]$ concentration of species $S$, $D \in \mathbb{R}_{+}$ diffusion coefficient, $v([S](t,x))$ convective velocity, and $v_{p}([S](t,x)), v_{c}([S](t,x))$ production and consumption rates.

Differential equations (especially ODEs) are the most commonly employed techniques to model biochemical systems due to their strong establishment in the sciences. Nevertheless using these methods (particularly PDEs) may also represent a significant mathematical challenge when attempting to solve large systems of non-linear differential equations. Moreover, it has been argued that the main challenge of this approach is the limited ability to describe biochemical systems with low species concentrations \citep{fontanaBarrier}. Chemical kinetic models specify the cell with limited structural descriptions. Biological systems are made of collections of objects whose identities are maintained and continuously evolve. These evolving properties may include the activation state, concentration, or the location.

\subsection{Markov chains}
Another method to examine biochemical systems is to express them as Markov chains \citep{gomez2001ppu}, in which
the state of the chain represents either approximations or exact number of the molecules present. Reactions are modelled as transitions between these states. The system is memoryless (``Markovian'') since the future development only depends on the present, not on the past. Therefore, the term Markov chain denotes time-discrete systems which are defined as a sequence of random variables $X_{1}, X_{2}, X_{3},...$ with the Markov property, i.e., 
$ P(X_{t+1} = x|X_{t}=x_{t},X_{t-1}=x_{t-1},...,X_{1}=x_{1}) = P(X_{t+1}=x|X_{t}=x_{t})$.

Provided there is no feedback in the system, the analysis of Markov chains is well developed, and the steady-state probability distribution of the process can be derived. Feedback, which is an inherent feature of many reaction networks, poses problems for analysis since a steady-state distribution of the system does not have to exist in this case.

Many straightforward, yet interesting simulation techniques which utilise the Markov property are based on explicit collisions between randomly selected molecules. This technique has the advantage of being easy to implement in a non-spatial case, and yet simple to extend to spatial simulations. A representative example of this type of algorithm is given by StochSim \citep{lenovere-shimizu-01}.

\subsection{Chemical master equation}
Where the model's time is continuous rather than discrete, the Markov chain is replaced by a ``continuous-time Markov process''. Here, the system again has a finite, discrete set of states, but now  a continuous time index $t$ exists. For simplicity, we focus on the case in which each state is given by the number of molecules per molecular species (i.e., a vector $x \in \mathbb{N}^k$). At any given point in time, the system occupies each state with a certain probability, yielding a probability distribution over all the states. The Chemical Master Equation (CME) provides a means to describe the temporal change of this distribution exactly for the case of a well-stirred and homogeneous reactor space \citep{vankampen2007spp}. Since chemical systems can be considered as Markovian, the CME approach is a special case of the continuous-time Markov chains. 

\cite{gillespie1976gmn,gillespie1977ess} proposed two precise ``Stochastic Simulation Algorithms'' (SSA) to simulate instances of the random process defined by the CME. These algorithms are widely used in the stochastic simulation of biochemical reactions \citep{meng2004mas} due to their significant efficiency in terms of computational cost. The principal factors in SSAs are reaction propensities $ f_{\mu} $, i.e., the likelihood of a reaction $ \mu $ to occur in the next (small) time step $dt$. These are computed from the mesoscopic rate constants and the number of molecules available as substrates to the reaction. From these, the next reaction and the time for that reaction have to be decided. This is done by using two random numbers. From the CME, it can be shown that the probability density function for reaction $ \mu $ to occur as the next reaction after time $ \tau $ is $ P(\mu,\tau) = f_{\mu} exp(-\tau\sum_{j} f_{j}), $ which is the basic equation SSAs are built on.

Gillespie's original work has been extended several times, most notably by the ``Next Reaction Method'' \citep{gibson-bruck-00}. This reduces the complexity from linear to logarithmic time in the number of reactions. Another technique is given by the ``tau-leap methods'' \citep{gillespie-01,chatterjee-vlachos-05}, which approximates the exact solutions obtained from SSAs. For larger numbers of molecules and reactions, however, these algorithms still suffer from high computational requirements. \cite{bernstein-05} extended the Gillespie algorithm to reaction-diffusion equations by dividing the reaction volume into several compartments and modelling diffusion between them.

\subsection{Bayesian networks}
A Bayesian network (BN) is a directed acyclic graph commonly used as a probabilistic modelling tool \citep{534975}. Modelling chemical networks with BNs was introduced by \cite{sachs02}. In a BN, variables (a molecular property) are represented as nodes in the graph. Directed edges express the dependence relation between nodes. A variable can be either discrete or continuous and may form a hypothesis, a known value (e.g., a concentration) obtained by experimental measurement or a latent variable. Variables which are not connected by edges are ``conditionally independent''.

If the state of a variable is known then the state of other variables can be predicted. This is accomplished through the use of: 

\begin{equation}
p(x)=\sum\ _y p(x,y)
\end{equation}
This formula sums the probabilities of all routes through the graph, thus allowing one to predict, with some probability distributions, the state of an unknown variable $x$. Continuous values for probabilities could be specified with a probability density function (e.g., \citealp{citeulike:461395} employs Gaussian distributions).

BNs have been used to reverse-engineer and infer the structure of biochemical networks \citep{sachs02,kim2003ign,citeulike:461395}. However, the setting of probabilities (learning) of BNs requires static experimental data, otherwise this may result in increasing the complexity of the task \citep{DBLP:journals/jcb/LiL05,Ch96}. The solid foundation of BNs in statistics enables the handling of the stochastic behaviour of real chemical networks and noisy experimental measurements \citep{deJong2002mas}. Another attribute of using BNs is that they can be employed when incomplete or only steady-state data on the reaction network are available. In this common case, kinetic models have been found to be less useful \citep{citeulike:126817}. \cite{citeulike:470218} discussed the various techniques to infer BN models from experimental data.

\subsection{Term rewriting systems}
\label{sectionBLreview}
Regulated term rewriting is a basic principle of information processing. Biomolecules, their polymeric subunits or groups of similar biomolecules are interpreted as objects encoded by character strings (terms). Sets of term rewriting rules describe possible interactions among objects and system components (e.g., pathways or membrane structures). Each application of a rule performs a discrete step of a process. The terms as a whole contain all information about the system status. Term rewriting systems can run in a massively parallel manner considering nondeterministic recombinations. Classes of grammar systems, P-systems \citep{paun1998dcn}, broadcast language \citep{holland1975ana,holland1992ana} and Alchemy based on the lambda calculus fall into this category \citep{fontana1994aft}. We demonstrate this modelling approach with the broadcast language (BL).

Holland originally proposed the BL formalism to assist his research on the ``adaptive plan''. Holland argued that the BL provides a straightforward representation for a variety of natural models such as biochemical networks.

The BL basic components are called \emph{broadcast units} which are strings formed from the set of ``monomers'' $\Lambda=\{0,\ 1,\ *,\ :,\ \lozenge,\ \triangledown,\ \blacktriangledown,\ \triangle,\ p,\ '\}$. Molecular species are broadcast units which can be viewed as \emph{condition/action} rules. Whenever a broadcast unit conditional statement (pattern matching expression) is satisfied, the computational action statement is executed, i.e., when an enzyme broadcast unit detects, in the environment, the presence of one or more specific substrate signal(s) then the broadcast unit broadcasts an output product signal. General signal processing can also be performed with broadcast units: e.g., a broadcast unit may detect a signal $I$ and broadcast a signal $I'$, so that $I'$ is some modification of the signal $I$. The broadcast monomers/symbols encode for the pattern matching and computational/enzymatic functions of molecular species. In addition, broadcast symbols may act as both operators and operands addressing the reflexive nature of molecular species (i.e., a molecule may act as both an enzyme and/or substrate).

Limited stochastic elements are involved in the computational functions of broadcast units which result in a semi-stochastic behaviour at the system level. The modelling of a genetic regulatory networks (which addressed only the regulatory/qualitative aspects of CIPNs) using the BL was proposed by \cite{DBLP:conf/eurogp/DecraeneMMK07}. Although possible, no quantitative studies have been previously reported to have been conducted with the BL prior to the work described in this thesis. Finally the BL formalism does not account for spatial information.
 
\subsection{Petri nets}
Petri nets (PNs) are a graph-oriented formalism originally from formal software engineering. Developed in the early 1960s \citep{Petri62,539513}, Petri nets provide a means to model and analyse systems, which comprise of properties such as concurrency and synchronisation. Petri nets consist of ``places'', ``transitions'', and ``arcs''. ``Arcs'' are used to connect the ``transitions'' and ``places'', ``input arcs'' connect ``places'' with ``transitions'', while ``output arcs'' start at a ``transition'' and end at a ``place''.

The modelling of biochemical networks with Petri nets was introduced by \cite{reddy1993pnr}. Here, place nodes are used to represent molecular species (enzymes, compounds, ions etc.) and transition nodes to denote chemical reactions. Other elements can be defined to specify in detail the chemical reactions to occur \citep{citeulike:273}.

Ordinary Petri nets provide an accessible modelling tool with well-established analysis techniques. For this reason, the use of Petri nets for qualitative analysis of biochemical network is growing. However, due to their timeless nature, Petri nets are limited regarding dynamic network analysis.

\subsection{\texorpdfstring{$\pi$}{pi}-calculus}

\newcommand{\emptyprocess}{0}

The $\pi$-Calculus is a process calculus, which is a formal method for modelling concurrent communicating processes \citep{358021,329902}. The $\pi$-Calculus provides a framework for the representation, simulation, analysis and verification of such systems.  The $\pi$-calculus allows the application of algebraic reasoning in order to determine the  equivalence between processes.

When modelling biochemical networks using $\pi$-Calculus, molecules and their individual domains are treated as computational concurrent processes \citep{DBLP:conf/psb/RegevSS01}. Complementary structural and chemical determinants correspond to communication channels. Chemical interactions and
subsequent modifications coincide with communication and channel transmission.  

The $\pi$-Calculus provides a highly detailed description of network nodes. However, the basic $\pi$-Calculus gives only a semi-quantitative view. A  significant factor to be considered is the lack of an associated temporal dimension and as a result all interactions can occur with the same
probability/rate. Extensions of the basic $\pi$-calculus address this limitation \citep{regev2004cab,blossey:17}.

\subsection{Agent-based models}
\label{subsec:ABLCS}
In an agent-based model (ABM), several computational objects called \textit{agents} are simulated to reproduce real phenomena within an artificial environment. ABMs originate from the late forties with the 
development of Cellular Automata \citep{vonneumann1949tao} and have been extensively used in the following fields: complex systems, multi-agent systems, and evolutionary programming \citep{994670, 975254}. An ABM is typically implemented with an object-oriented framework \citep{rumbaugh1991oom,bersini2004ooo}. Each agent or class is defined with particular properties and methods. Agents are situated in space and time, interactions between with each other may occur following rules. Global and complex behaviour may emerge from these local agent-agent interactions and properties. 

ABMs provide a flexible framework to: specify and refine with ease rules governing agent behaviours and interactions (e.g., using production rules or Boolean logic), secondly, to model emergent system or global behaviours \citep{pmid16223502}. Preliminary works to model bio-chemical networks using ABMs appeared in the late nineties \citep{pmid9076693, pmid10400267}. ABMs consider the cell and its components as agents with cognitive capabilities. 
Two distinct ABM approaches are presented:
\begin{enumerate}
	\item In Cellulat, which was developed by \cite{Gonzalez2002mis,gonzalez2003cab}, a cell is seen as a collection of adaptive autonomous agents. Communication between agents is performed via propagating signals on a shared data structure, named ``blackboard'' referring to the blackboard architecture \citep{DBLP:journals/aim/Nii86, DBLP:journals/aim/Nii86a}. An agent receives a signal or a combination of signals from a designated blackboard level and transduces these into another signal (or set of signals) on the same or different blackboard level. Transduction mechanisms of the signal depend of the cognitive capabilities of the agent. A blackboard level could represent extracellular, membrane, cytosol or nucleus region, this enables the modelling of spatial organisation.

	\item A second ABM is described where Learning Classifier Systems (LCS) are used to specify the agents' behaviour and interactions. LCS are systems constructed from condition-action rules called \textit{classifiers}. LCS can be seen as a simplification of the broadcast language where classifiers are binary strings that can be viewed as IF/THEN statements. Holland's initial work was modified a number of times and at present many different varieties of learning classifier systems are available \citep{2001-Lanzi-0, bull2005flc}. 

LCS are commonly used as a machine learning technique. However \cite{holland2001eec} proposed an agent-based model where the agents' behaviour and adaptation are determined by the use of LCS. This work argued that LCS could be used to evolve a simple repertoire of condition-action rules to a more complex goal directed set of rules. 

In typical biochemical networks, interactions between molecules follow the same condition-action mechanisms. Thus Holland suggested that this approach could be used to model and simulate CIPNs. His proposition to design chemical networks was to start with a LCS-based ``over-general" model of a biological phenomenon (e.g., transformation of a healthy cell to a cancer cell). Then this general phenomenon could be refined through several iterations. At each iteration, the details (e.g., compartment level) of the occurring interactions can be specified. These iterations were continued until the desired network level/granularity was reached, where the submolecular objects are specified (e.g., protein ligand, receptor, ions etc.). This refining process highlights the top-down/hierarchical approach and descriptive power of LCS to model and simulate complex CIPNs. Moreover this approach can be naturally coupled with Genetic Algorithms. This evolutionary feature may allow one to examine phylogenetic relationships between different reaction networks (where the signalling differences may be due to random molecular mutations).
However no actual implementation and experimental examination of this system have ever been reported, therefore this proposal and associated potential benefits remain conjectural.
 
\end{enumerate}

\subsection{SBML \& CellML}
Modelling techniques may be employed in conjunction with a markup language to store generated models. The use of a standard format facilitates the analysis, visualisation, simulation and exchange of biochemical network models within the modelling community, providing opportunities for refinement and incorporation of new knowledge. So far, two approaches have emerged, resulting in the model-description languages SBML (Systems Biology Markup Language) \citep{hucka2004elf} and CellML \citep{lloyd2004cif}, both based on the XML markup language \citep{bray2000eml}.

\begin{itemize}
	\item In SBML, a biochemical network is described in terms of the molecules taking part in it - termed species - and the reactions taking place between them. The present amount of each species can be expressed either in terms of its concentration or of the number of molecules present. Each reaction has an associated kinetic law, which defines the rate of the reaction depending on the present amount of its substrates. Additionally, the model can be subdivided into a fixed set of well-stirred compartments to include a non-hierarchical spatial component. Nevertheless SMBL models cannot specify fluxes between compartments at present (i.e., in SBML level 2 version 4 release 1).
	\item  In CellML, a more general approach is taken, in which a model consists of components and connections between components. Each component can contain variables and a reaction between them, and connections are used to transfer the value of variables from one component to another.
\end{itemize}

Although CellML is following a slightly more general approach, it is not as widely used as SBML, for which a large collection of software tools is available (see www.sbml.org for a list of these tools).
Additionally, the first model repositories have started to use SBML as a representation language, e.g., see the BIOMODELS database at \url{www.ebi.ac.uk/biomodels}. Therefore, SBML can be seen as the first emerging specification standard for biological models at the cellular level. Finally the use of such a common language provides the ability to analyse and complement intersecting information on differing compatible modelling techniques. 

\section{Comparison of approaches}
In this section, we compare the previously introduced methods to model CIPNs by using a set of defined criteria. Following this, a comparison table is presented to summarise this review. The intention is to determine a suitable modelling technique to be employed in our project. This selection will be discussed after the presentation of this comparison table. We identify evaluation criteria with regards to stochasticity, time, granularity, space, topology and modularity.  

\subsection{Evaluation criteria}
Relevant criteria are outlined here in order to compare the modelling techniques presented in Section \ref{sec:surveyModel}:
\begin{itemize}
	\item \textbf{Stochasticity}: This property reflects the range of possible processing scenarios that may be identified by the model.
  			\begin{itemize}
					\item \textit{Deterministic}: The system behaviour purely depends on inherent data. No external or statistical fluctuation may occur and influence the system's dynamics. The system may only operate along one known path.
     			\item \textit{Nondeterministic}: A number of alternative paths for system processing may exist which can be completely explored by the model. All possible scenarios are taken into account by the model in which no unanticipated events may affect the system's dynamics. 
     			\item \textit{Stochastic}: In contrast, stochastic models select one possible path in a random manner that can be based on a given probability distribution. This implies uncertainty (external and statistical fluctuation may be accounted for) and inhibits repeatability of systems runs.
  			\end{itemize}
	
	\item \textbf{Time}: This property describes how time is represented within the model.
				\begin{itemize}
					\item \textit{Atemporal}: When executed, the model remains static and introduces no temporal consideration.
					\item \textit{Events}: A sequence of pre-identified events defines the granularity of time. An event is an action within the system which characterises the progress of the system processing. Events are not necessarily equidistant in time. Dependencies between processes, their synchronization and concurrency may also be based on the interplay of events.
					\item \textit{Discrete}: Temporal changes are characterized by fixed periodic intervals. A discrete time interval defines the smallest unit measuring the system's dynamic behaviour. Discrete time points allow one to express recursive formulation of the system processing. Discrete time may be referred as a global clock for the system.
					\item \textit{Continuous}: Infinitesimal time intervals allow the finest granularity for measuring time represented by real numbers. Computer-based simulation techniques, by their nature, require an approximate discretisation of points in time.		
				\end{itemize}
   
   
   \item \textbf{Granularity}: This property designates how the molecules or particles are represented in the model. It refers to the abstraction level of their specification. The finer the granularity the more detailed the system that can be described. Granularity also constrains the level of monitoring capabilities.
 
				\begin{itemize}
					\item \textit{Submolecular}: This level allows one to compose molecules by atomic specifiers or functional units (e.g., protein domains).
					\item \textit{Molecular}: Molecules are considered as the smallest expressible unit. A mapping between the chemical substance and the assigned identifier (e.g., symbol) is either assumed or abstracted.
					\item \textit{Species}: An enumerable amount of molecules having the same chemical substance is regarded as a species. This level of granularity enables one to quantify a molecular species as a whole within the system, however one cannot isolate an individual molecule of a given species.
					\item \textit{Concentration}: Allows one to quantify the relative amount of a particular molecular species existing in a system. As represented by real numbers, transforming absolute molecular amounts into concentrations can require an approximation. Concentrations can be viewed as an approximation of the molecular species quantities.
				\end{itemize}
   

	\item \textbf{Space}: When handling molecules of given granularity within a model, a system component which is analogous to a reactor is assumed. This component can provide space if the positioning of the molecules (within the reaction system) is taken into consideration. 
				
				\begin{itemize}
					\item \textit{Implicit}: Particle or molecule identifiers include spatial information, e.g., using an index. System components that control the evolution can be equipped with regulation schemes for updating this information. Here, a homogeneous distribution of the molecules within the reactor is assumed. In this ``well-stirred'' reactor, no boundaries are specified, and there is no explicit definition of space in the model.
					\item \textit{Compartmental}: A hierarchically nested or graph-based number of explicit compartments is distinguished. Each molecule is assigned to one of the specified compartments and can move from one compartment to another. Within each compartment, no further specification of molecular positioning is defined.
					\item \textit{Grid}: Apart from the compartmental structure, a spatial geometry is used to locate molecules more precisely. This way, discrete spatial distributions of molecules can be mapped using the model.
					\item \textit{Continuous}: The finest granularity of defining space is given by continuous values. Here, each molecule can be positioned arbitrarily within the reactor. Analogous to continuous time, computer-based simulations may require discretisation which would imply approximation.
				\end{itemize}


	\item \textbf{Topology}: This designates the ability of the model to dynamically modify its structural components (e.g., pathway structure, dependencies between compartments, active membranes, receptor dynamics). 
 
				\begin{itemize}
					\item \textit{Fixed}: A static system structure is assumed.
					\item \textit{Dynamic}: Principles or rules are defined that allow the system structure to change over time and space. These rules are a part of the model description.
				\end{itemize}
  
	\item \textbf{Modularity}: This refers to the ability of the model to subdivide a given biological reaction system into functional sub-units (i.e., modules). The subdivision process is carried out through algorithmic strategies applied on the model. Modules are determined/classified according to specific properties (e.g., network topology/clusters, functions) across these sub-units. Modularity may facilitate the study of a system by examining sub-units independently instead of the system as a whole.
				\begin{itemize}
					\item \textit{No}: The whole reaction system is regarded as a monolithic entity which currently prevents the identification of sub-units.
					\item \textit{Hierarchical structure}: The sub-units are represented as nodes forming a tree-based structure. Modules communicate with each others (e.g., transmission of molecules from one sub-unit to another) via specified interfaces, typically through diffusion over transduction/communication channels.
					\item \textit{Graph-based structure}: These structures are a generalisation of tree-based structures which does not necessarily account for a hierarchical organisation.
				\end{itemize}
\end{itemize}

\subsection{Comparison table and discussion}
As a summary of previous sections, a comparison table is presented (Table \ref{tab:revModel}) which uses the criteria that were discussed above. The table provides an immediate comparison of differing modelling techniques and allows one to identify desirable attributes which may be necessary for modelling a specific biochemical system.

\begin{table}[htbp]
\centering
\rotatebox{90}{
{\scriptsize
\begin{tabular}{|l@{}||c@{}|c@{}||c@{}|c@{}||c@{}||c@{}|c@{}|c@{}||c@{}|c@{}|}
\hline
 & \rotatebox{90}{Differential equations (ODE)} & 
 \rotatebox{90}{Differential equations (PDE)} & 
 \rotatebox{90}{Markov chains} & 
 \rotatebox{90}{Master equation (CME)} & 
 \rotatebox{90}{Bayesian networks} & 
 \rotatebox{90}{Term Rewriting Systems (BL)} & 
 \rotatebox{90}{Petri nets} & 
 \rotatebox{90}{$\pi$-calculus} & 
 \rotatebox{90}{Cellulat} &
 \rotatebox{90}{Agent-based LCS}
 \\
\hline
\hline
Model type & deterministic & deterministic & stochastic & stochastic & probabilistic & algebraic & algebraic & algebraic & agent-based & agent-based
\\
\hline
Time & cont. & cont. &  discrete & cont. &  atemporal&  discrete & events  & events & discrete & discrete \\
\hline
Granularity & conc. & conc. & species. & species & species & (sub)molecular & molecular  & molecular & (sub)molecular & (sub)molecular  \\
\hline
Space & implicit & cont. & implicit &  implicit&  implicit& implicit & implicit  & implicit & compartmental & implicit\\
\hline
Topology & fixed & fixed & fixed & fixed & fixed &  dynamic& fixed  & fixed  & dynamic & dynamic\\
\hline
Modularity & no & no & graph-b. & graph-b. & graph-b. &graph-b. & graph-b.  & graph-b. & $\begin{array}{l}\textrm{hierarch. or}\\\textrm{graph-based}\\\end{array}$ & hierarch. \\
\hline
\end{tabular}
}
}
\caption{Comparison of CIPN modelling approaches with respect to previously defined classification scheme}
\label{tab:revModel}
\end{table}
 
\clearpage
Based on our review of modelling techniques which is summarised in Table \ref{tab:revModel}, we discuss and relate this evaluation with the requirements of this project identified as follows:

\begin{itemize}
	\item The specification and development of a software platform capable of modelling and simulating artificial CIPNs.	
	\item The examination/traceability of the individual behaviour of molecular species/instances.
	\item The possibility to conduct both qualitative and quantitative analysis of the CIPN system's dynamics.
	\item The inclusion of stochastic processes reflecting the random nature of molecular collisions/reactions.
	\item The specification of submolecular properties (e.g., monomers).
	\item The ability to dynamically change the CIPN's topology over time (e.g., creation/modification or removal of molecular reactions due to molecular mutations).
\end{itemize}

Probabilistic models are employed to infer the model from experimental observations which conflicts with our attempt to construct a simulation platform. Examining organisational closure in CIPNs implies that the monitoring/tracing of the individual molecules is required. However deterministic and stochastic techniques treat the molecular species in an aggregate manner. As a result algebraic and agent-based approaches are more appropriate for the current endeavour. 

On the other hand, the main drawback of algebraic and agent-based approaches is the fact that they may not always permit a detailed (i.e., both qualitatively and quantitatively) analysis of biochemical processes' dynamics. The latter is a key feature of deterministic and stochastic approaches. Moreover, these mathematically grounded approaches have a strong theoretical scientific foundation associated with a plethora of analytical and simulation tools available. 

Nevertheless we may also consider that we intend to realise and evolve \textit{artificial} CIPNs. In other words, we do not intend to provide a simulation platform of \textit{real} biochemical networks, in which case a precise account for the temporal dimension would have been necessary. Thus the semi-quantitative and discrete approach of algebraic and agent-based techniques is sufficient. Moreover, exporting algebraic/agent-based models into the SBML format allows us to conduct further analyses of the systems' dynamics using complementary SBML tools. In Chapter \ref{chap:evolClosReacNet}, we employ the SBML to generate and analyse the deterministic dynamics of reaction networks employed and generated in our evolutionary experiments.

Finally agent-based approaches offer most flexibility: ABMs provide a detailed and adaptable description of molecular species (including submolecular components) and of the system as a whole. Secondly it is possible to utilise existing algebraic methods (and exploit associated features) to specify molecular species and reactions within an ABM. Most ABMs account for the dynamic nature of biochemical systems in which the relationships between molecular species/modules may change over time due to internal or external perturbations (e.g., molecular mutations). Stochastic elements may also be involved in the specification of the agents' interactions with other species/systems' elements.

According to the requirements of this project, ABMs appear to represent a suitable and flexible technique to model, simulate and analyse CIPNs and their evolution. 

\section{Conclusion}
We introduced the overall concepts associated with the main CIPN modelling branches: deterministic, stochastic, probabilistic, algebraic and agent-based approaches. To best illustrate each of these classes of modelling techniques, we presented and reviewed a selection of techniques.
 
We also identified two markup languages: SBML and CellML which are widely used to specify, disseminate and exchange CIPN models within the scientific community. These languages also permit the migration of specifications between differing modelling techniques. This migration capability allows one to widen the range of analytical studies of a given modelling technique. We noted that the SBML possesses a longer history than CellML and has subsequently become the standard language for storing CIPN models. We will therefore employ the SBML as a means to migrate and disseminate our generated CIPN models.

Following this, a model comparison table was presented and highlighted the modelling capabilities of the different techniques. We outlined the requirements of this project with regards to the modelling of CIPNs. We related and discussed these requirements with the review conducted in this chapter. We finally distinguished ABM approaches as a suitable and flexible modelling technique. We consequently select the ABM framework to model, simulate and analyse CIPNs in this thesis.
