...

SIMULATION OF A HETEROGENEOUS BOOLEAN NETWORK

by user

on
Category: Documents
12

views

Report

Comments

Transcript

SIMULATION OF A HETEROGENEOUS BOOLEAN NETWORK
1
SIMULATION OF A HETEROGENEOUS BOOLEAN NETWORK
Celeste Mott
University of Nebraska at Omaha
Mathematics
2012-2013 KRMP and UNO FIRE grant
Advisor: Dora Matache
ABSTRACT. A Boolean network is a system of nodes linked to one another, each
node capable of two states, active or inactive. As time passes, each node examines
the state of its input nodes and determines its own state according to a specific rule,
or Boolean function. As the network evolves, it can either regress into a cycle of
repeating states or erupt into chaos under certain conditions. This paper examines a
computer simulation that creates a Boolean network of n nodes where each node
has a connectivity (number of input nodes), k, drawn from a uniform distribution,
and one of four general function types—canalizing with one canalizing input,
canalizing with two canalizing inputs, threshold, or bias functions. The goal is to
estimate the probability that the network returns to its original behavior after a
small disturbance. The main purpose of this computer simulation is a means of
reference. Since the program has a great number of editable parameters, one can
use the simulation to estimate how modifying a parameter will affect the stability of
a network or compare that with mathematical models or real data against the
simulation.
1. INTRODUCTION
A Boolean network is a collection of nodes that have two possible states, active or
inactive, alive or dead, 1 or 0. Each node in the network is connected to other nodes in the
network and / or itself. The node examines the state of its input nodes, then activates or
deactivates at the next time step according to a specific rule, or Boolean function. Many
biological networks can be described as a Boolean network. For example, in a cellular
network, protein A may activate when in contact with an active protein B or it may
deactivate when in contact with an active protein C.
The main goal when studying a Boolean network is to understand the long term
dynamics of the system and identify if the network is stable or not. A stable network
reaches a cycle of states and an unstable network evolves in an unordered fashion known
as chaos, which is caused by inherent disturbances.
Most of the work in Boolean network literature has examined networks for which
all nodes obey the same Boolean function. However in real networks, such as genetic signal
2
transduction networks, the network is composed of nodes with several different Boolean
functions. The proportion of nodes following each type of Boolean function varies and
nodes from one category may prefer to be linked of nodes of a different category.
For example, in [Derrida_Pomeau] the authors identified the so-called critical
connectivity value, k, which separates order from chaos under a variety of network
parameters for a biased Boolean network. Biased Boolean functions generate an output of
1 with a given probability, ρ, called the bias. The critical condition is given by the following
equation which describes the average influence or sensitivity of the network, or the
probability that after changing the state of a random unit of the network in a limiting cycle
the system returns to the same cycle.
(
)
(1)
Figure 1 was obtained in [Aldana_Cluzel] and shows how the parameter space k-ρ is
divided into two regions—order (grey) and chaos (white). The curve separating the two is
the critical curve given by equation (1).
Figure 1. (Caption from [Aldana_Cluzel]) Dynamical robustness of a Boolean
network with random topology. The network dynamics exhibit both chaotic and
robust behaviors depending on the value of the parameters k and ρ. The parameters
space k-ρ is then divided into two distinct regimes—chaotic (white) and robust
(grey). The curve separating these two regimes is given by the solutions of (1). The
parameters space is total dominated by the chaotic regime.
3
Similarly, in [Moreira_Amaral] the authors find the phase diagram which indicates the
phase transition from order to chaos and the critical curve for a network governed by
canalizing Boolean functions with one canalizing input. A Boolean function is canalizing if,
whenever one variable (the canalizing variable) takes a given value (the canalizing value)
the function always yields the same output (canalized value). As indicated in [MoreiraAmaral], studies of gene regulation in eukaryots have showed that the Boolean idealization
is a good approximation for the nonlinear dynamics of this system and that the gene
regulating mechanism have a strong bias toward canalizing functions [Harris_Etal]. The
critical condition that separates order from chaos is obtained by computing the average
influence/sensitivity of the network and is given by the following equation from
[Moreira_Amaral]
(
)
(
)
(
)[
(
)] (
)
(2)
where ρ1 = fraction of functions that have canalized output 1
ρ2 = bias of the function that takes over when the canalizing input is not on its
canalizing value
η = fraction of functions with canalizing input 1 and q is the long-run average
fraction of nodes in the active state, which can be estimated from the
known parameters as the solution of the system
(
where
(
)(
)
)
P = probability that the canalizing input is on its canalizing value
Figure 2 is obtained from [Moreira_Amaral] and represents several phase transition
diagrams corresponding to networks of increasing size N and with variation of parameters
as specified in the graphs.
4
Figure 2. (Caption from [Moreira_Amaral]) Order-chaos phase transition for
networks of canalizing Boolean functions. The color coding represents the
probability that after changing the state of a random unit of a network in a limiting
cycle the system returns to the same cycle. The resilience of the network to small
perturbations is a signature of the ordered phase, corresponding to small
probabilities. The chaotic region is mainly dark blue and the critical curve given by
equation (2) is plotted with a dotted black curve.
Recent work [Ludes_Mott_Matache] has been done to identify critical conditions for
networks governed by canalizing functions with two canalizing inputs, such as whenever
two inputs (the canalizing variables) take a given combination of canalizing values (00, 01,
10 or 11) the function always yields the same canalized value (0 or 1). The critical
condition is given by
[ (
(
where
(
(
(
)(
)
(
)[ (
)[(
)[
)[(
)]
)
)(
(
)
(
)
)
)
(
(
)]
)]
(
(
)]
)]
(3)
q = the long-run average fraction of nodes in the active state
= ( ∑
{
)
}
(
)
5
P = probability that any of the canalizing inputs is on its canalizing value
(
=
)
(
)
)(
k = connectivity
ρ = bias of the function that takes over when the canalizing input is not on
its canalizing value
= the proportion of nodes in the network with canalizing input values xy
and canalized output z where
{ }
Similarly, in [Ludes_Mott_Matache] the critical condition is obtained for a network
governed by Boolean threshold functions such that a node becomes active if and only if the
fraction of active inputs with value 1 is within some thresholds d1 and d2 with 0 ≤ d1 ≤ d2 ≤
1; otherwise the node is deactivated (given state 0). The critical condition is
[∑
(
(
)
) [∑ (
(
)
∑ (
)) ]
∑ (
(
)) ]
(4)
where q = long-run average faction of nodes in the active state obtained from the map
(
)
( ) [∑ ( ( )
) ( )
(
where
α =
=
γ =
=
(
(
(
(
)∑
∑ ( ( )
( ( )
)(
( ))]
)
N = number of nodes in the network
(
) ( ( )) (
)
is a binomial term
|(
|(
|(
|(
|
|
|
|
)
)
)
)
)
)
)
)
where σ1 = canalizing input that is flipped from 1 to 0 or from 0 to 1 as
specified
The quantities α and γ are obtained by considering all possible combinations of d1 and d2.
The work in [Ludes_Mott_Matache] is in the process of being supplemented with phase
transition diagrams using the work in this report, which examines a computer program
6
written to generate and simulate Boolean networks in order to determine the phase
transition diagrams and critical curves for networks governed by canalizing functions with
one canalizing input and networks governed by threshold functions. A further goal of this
work is to provide a tool for simulating the phase diagram for a heterogeneous network
which combines all different types of Boolean functions mentioned above, including biased,
canalizing with one or two canalizing inputs and threshold. As a matter of fact, the scheme
shown in Figure 3 describes the overall goal related to biological networks emphasizing
with colored boxes the focus of this particular report. The overall work involves both
providing a mathematical model for the critical condition which is explored in detail in
[Ludes_Mott_Matache] and comparison with the results of a simulated network with
similar characteristics, which is the goal of this report.
Figure 3. Scheme of the heterogeneous Boolean network model; the mathematical
model is used to generate the critical condition and matched with the phase diagram of
a simulated network. The focus of this report is highlighted in pink boxes and is related
to the simulated network which is to be compared to the prediction of the model.
7
In the next sections we present the computer algorithm used to simulate the
heterogeneous network. All work was written in the language C++ and GUI development
used the [Qt] library. The complete code is presented in the appendix.
2. INITIAL NETWORK SETUP
The program begins by launching a GUI dialog. Figure 4 shows that the dialog contains a
number of parameters that are adjustable by the user, beginning with number of nodes, n.
The number of nodes can be manually selected or randomly generated. Currently the
maximum number of nodes in the network is 52 because of the way the program identifies
network cycles but the maximum number can be expanded with some work.
Next the user may select connectivity k. Like n, k may be manually selected or randomly
chosen between 1 and 5. The user may also select if each node is to have exactly k inputs or
between 1 and 2k-1 inputs for a network average connectivity of k.
Finally the user specifies the parameters for each function. There are four possible
types of functions including threshold function, canalizing with one canalizing input,
canalizing with two canalizing inputs or bias function. The user must select what
percentage of nodes in the network fall under each category. The program is safe-guarded
so that the network will not generate if the sum of the percentage of each function is not
100 percent or if the upper threshold of a threshold function is less than the lower
threshold. After that, specific parameters of each function type may also be adjusted—
please note these parameters are further explained in Simulating the Network, where we
describe how each function evaluates a node’s inputs and determines whether the node
should activate or deactivate.
After all the user’s parameters have been accepted, the program starts generating a
network by creating a list of nodes, n0, n1,… ,nn-1, that compose the Boolean network. Once
the list is created, each node is assigned a series of parameters beginning with function
type. The rules of each of the four functions are thoroughly explained under Simulating the
Network. The user specifies what percentage of nodes in the network should fall under
each of the four categories and the function of each node is randomly assigned to meet
these chosen proportions.
After the function type is assigned, the connectivity (number of input nodes), k, is
determined. The program currently has two options for determining individual
connectivity, a) all nodes have exactly k inputs or b) all nodes have approximately k inputs
for an average connectivity of k (using a uniform distribution—each node has an equal
chance of receiving 1 to 2k-1 inputs). Desired connectivity and distribution are parameters
both provided by the user.
After each node is assigned a connectivity, k, specific nodes are selected to be inputs.
Recall that each node is provided an index, i. The program randomly picks k numbers
between 0 and n-1 inclusive, then assigns the nodes with those indexes as inputs. Currently
the program does not allow single loopback, meaning if a node’s own index is randomly
generated that number is thrown out and a new index is randomly selected, however two
nodes may input to each other.
Once each node in the list is assigned a function, connectivity, k, and a list of inputting
nodes it is randomly assigned a state, 0 or 1. Currently any node has a 50 percent chance of
being active or inactive.
8
Figure 4. Screen shot of GUI interface that allows user to configure the Boolean
network.
9
3. SIMULATING THE NETWORK
After the user has satisfactorily entered all the parameters a second tab shown in Figure
5, the Simulate Tab, appears. This tab shows a summary of the network along with a few
important parameters. First is a list of parameter sets to vary. When examining a Boolean
network, it is helpful to create a 3-D plot that investigates order and chaos of the system.
This graph varies two parameters then plots the average sensitivity of the network, or
probability that a small perturbation does not affect the cycle the network will reach, of the
network on the z-axis. Under the Parameter Sets to Vary Section, if a checkbox is checked
then the program will create a text file that contains a matrix storing network sensitivity
that may be plotted in [Matlab]. T describes the number of iterations the network is
allowed to evolve.
After any parameters are modified the Simulate button may be pressed to start the
simulation. Recall that each node in the network is assigned one of four specific functions,
including bias function, canalizing with one canalizing input, canalizing with two canalizing
inputs or threshold function. The network simulates by evaluating each nodes’ function at
time t then activating or deactivating each node at time t+1 accordingly. In what follows we
show how this is done for each individual class of functions.
a. Evaluation of Nodes with Function: Threshold
A threshold function is arbitrarily assigned two values—lower threshold, d1 = 0.30, and
upper threshold, d2 = 0.60 (note that these default values may be altered by the user). For
each node governed by this type of rule the threshold function calculates the percentage of
input nodes that are active, p, and then compares it to d1 and d2. The node will activate at
time t = t+1 if d1 ≤ p ≤ d2 and deactivate otherwise.
b. Evaluation of Nodes with Function: Canalizing with One Canalizing Input
A canalizing function with one canalizing input has two parts. First, one inputting node
is randomly selected to be the canalizing input, its current state defined as σ1. If σ1 is on its
canalizing value then node will activate at time t = t+1. Table 1 outlines four possible
functions a node with one canalizing input may follow. If σ1 is not on its canalizing value a
bias function takes over. A random number between 0 and 1, r, is randomly generated. If r ≥
ρ2 = 0.50 (note that this arbitrary bias may be altered by the user) then the node will
activate at time t = t+1. If the canalizing node is not on its canalizing input and
then
node is deactivated.
Function Name Summary
CAN_POS
CAN_NEG
NOT_CAN_POS
NOT_CAN_NEG
Table 1. Description of four functions that fall under the canalizing function with
one canalizing input category.
10
Figure 5. Screen shot of GUI that allows user to begin simulation.
11
c. Evaluation of Nodes with Function: Canalizing with Two Canalizing Inputs
Similar to a canalizing function with one canalizing input, first two inputting nodes are
randomly selected to be the canalizing inputs, their current states defined as σ1 and σ2. If σ1
and σ2 are both on their canalizing values than node will activate at time t = t+1. Table 2
outlines eight possible functions a node with two canalizing inputs may follow. If σ1 and σ1
are not on their canalizing values a bias function takes over. A random number between 0
and 1, r, is randomly generated. If r ≥ ρ2 = 0.50 (note that this arbitrary bias may be altered
by the user) then the node will activate at time t = t+1. If the either canalizing nodes are not
on their canalizing inputs and
then node is deactivated.
Function Name
CAN1_AND_CAN2_POS
CAN1_AND_CAN2_NEG
CAN1_AND_NOT_CAN2_POS
CAN1_AND_NOT_CAN2_NEG
NOT_CAN1_AND_CAN2_POS
NOT_CAN1_AND_CAN2_NEG
NOT_CAN1_AND_NOT_CAN2_POS
NOT_CAN1_AND_NOT_CAN2_NEG
Summary
Table 2. Description of eight functions that fall under the canalizing function with two
canalizing inputs category.
a. Evaluation of Nodes with Function: Bias
A bias function generates a random number, r, between 0 and 1, and compares it to an
arbitrary threshold, b = 0.50 (note that this default value may be altered by the user). The
node is activated at time t = t+1 if r ≥ b.
b. Averaging Results
In order to obtain the phase diagram we compute the average sensitivity, or probability
that after changing the state of a random node in a limiting cycle of states the system
returns to the same cycle of states, for the network. We estimate this probability using the
procedure described below. The critical condition is obtained when the average sensitivity
is set equal to 1/k for connectivity k. As shown in [Ludes_Mott_Matache] for the
heterogeneous network we use a weighted connectivity where the weights represent the
proportions of each class of Boolean functions in the network.
The network is evolved from time t = 0 to either time t = T (arbitrarily assigned by user)
or until a cycle is reached. A cycle occurs when the network reaches a state is was
previously in and as a result it will follow the same series of evaluations and continue in a
loop of states. (Note that evaluating to a cycle is not entirely accurate since the outcome of
a bias function is not dependent on a node’s current state or a node’s inputs.) Whether a
network reaches a cycle or not may depend on the initial state of each node, so to gain a
12
good approximation of the network’s behavior several simulations are run. First the
network is simulated until time t = T or a cycle is reached, then the state of each node i is
inverted (one at a time) to see if the same outcome occurs. This process takes 2n
simulations for a network with n nodes. However, flipping the state of each node may not
be enough to average out any anomalies so the process is repeated one hundred times for
one hundred different initial states of the network. Note that during each simulation, each
node retains the properties: function type, connectivity, k, and list of input nodes, the only
parameter that is being changed is the initial state—active or inactive.
After 200n simulations, a number, S, called the stability is produced. Stability, S, is
between 0 and 1 and describes the percentage of times that inverting a node’s initial state
had no effect on the network (the same cycle was reached as before, either a true cycle or
the collection of network states for t between 0 and t). Then 1-S is an estimation of the
probability that after changing the state of a random node in a limiting cycle of states, the
system returns to the same cycle of states. In other words it is an estimation of the average
sensitivity of the network whose values indicate the phase of the system: order, chaos, or
criticality (also known as edge-of-chaos) which produces the more complex behavior in the
system.
The process is repeated for a large number of parameter combinations and phase
diagrams can be obtained by plotting the value of 1-S against two selected parameters
which are varied over an entire interval. The critical curve is obtained for 1-S = 1/k, where
k is the weighted connectivity described previously.
The program generated provides a collection of output files that allows users to graph
the phase diagrams in a separate environment. In particular for the work in
[Ludes_Mott_Matache] the results are plotted with the technical computer software Matlab
[Matlab].
c. Analyzing Results
The main purpose of this computer is a means of reference. The goal is to examine how
parameters affect the stability of the network. Figure 6 shows an example of a text file
produced by the program that varies the parameters ρ1 and ρ2 of a canalizing function with
one canalizing input. The data is organized into a matrix such that the horizontal rows
represent different values for ρ1 between 0 and 1, the vertical columns represent different
values for ρ2 between 0 and 1 and the value represents the stability of the network. This
matrix of values can then be plotted by [Matlab] as shown in Figure 7.
13
Figure 6. Example of a text file produced by this program that displays the stability of a
network for 0 ≤ ρ1 ≤ 1 and 0 ≤ ρ2 ≤ 1.
Figure 7. Graph plotted in [Matlab] using text file produced by network simulation. The
graph clearly shows for which values of ρ1 and ρ2 the network remains stable or
becomes chaotic.
14
4. CONCLUSIONS AND FURTHER RESEARCH
This work provides a computer program for simulation of heterogeneous Boolean
networks and computes the average sensitivity of the network for a variety of parameters
in order to generate phase transition diagrams that indicate the transition from order to
chaos in terms of the underlying parameters. The program can be used to run a variety of
scenarios that may represent real network characteristics in order to determine which
parameters can be modified in order to drive the system to a desired behavior.
One immediate purpose of this work is to validate the associated model obtained in
[Ludes_Mott_Matache]. Once the model is validated it can be used independently of this
program for a variety of network scenarios. We plan on comparing the results obtained
with the simulator generated by this work and the model to the results obtained in
[Helikar] for a generic cellular fibroblast signal transduction network. In [Kochi_Matache]
the network [Helikar] is divided into several different classes of Boolean functions and it is
shown that the network exhibits order. The network consists of a fairly large number of
nodes connected by regulatory mechanisms. The Boolean rules of that network encompass
canalizing functions with one or more inputs (with complete or partial canalization; partial
canalization is subject for future research), threshold functions, while the remaining
functions can be regarded as biased functions. We will identify the parameters of that
network to assess the accuracy of our approach in estimating the behavior of that network.
Then we plan on altering different aspects of the network to study related networks that
may be of interest in genetics.
Although this work is targeted to biological networks, our procedure is valid for any
types of networks for which we can establish a partition of the ensemble of Boolean
functions. Thus, it is applicable to biological, chemical, physical, social, financial networks
etc.
The following is a list of ways to expand the program of this report:
1. Each node currently has a connectivity of exactly k or connectivity between 1 and
2k-1 for a network average of k. Distributions other than uniform could be applied
to the network.
2. In a biological network, certain proteins may have more or less connections than
other proteins. This program could be modified so that nodes of one function are
likely to have a different connectivity than nodes of a different function.
3. Currently each node’s input nodes are selected at random. This program could be
modified so that nodes of function type A are more likely to connect to nodes of
function type B than of type C, etc.
4. Currently each node has a 50 percent chance of being active at time t = 0. This
program could be modified so that nodes of function type A are more likely to be
initially active than nodes of function type B. Various initial distributions could be
explored.
5. More types of functions could be added to the network.
15
References
[Aldana_Cluzel] Aldana M., Cluzel P., A natural class of robust networks, PNAS 100(15),
2003, p. 8710-8714.
[Derrida_Pomeau] Derrida B., Pomeau Y., Random networks of automata: a simple
annealed approximation, Europhys. Lett. 1, 1986, p. 45-49.
[Harris_etal] S. E. Harris, B. K. Sawhill, A. Wuensche, and S. Kauffman, A Model of
Transcriptional Regulatory Networks Based on Biases in the Observed Regulation Rules,
Complexity 7, 2002, p. 23-40.
[Helikar] Helikar T., Konvalina J., Heidel J., Rogers J.A., Emergent decision-making in
biological signal transduction networks, PNAS 105, No. 6, 2008, p. 1913-1918.
[Kochi_Matache] Kochi N., Matache M.T., Mean-Field Boolean Network Model of a Signal
Transduction Network, Biosystems, Volume 108, Issues 1-3, April-June 2012, Pages 14-27.
[Ludes_Mott_Matache] Ludes A., Mott C., Matache D., Phase Diagram of Heterogenous
Kauffman Networks, preprint.
[Matlab] Matlab is a product of the Mathworks, http://www.mathworks.com/index.html
[Moreira_Amaral] Moreira A.A., Amaral L.A.N., Canalizing Kauffman networks:
nonergodicity and its effect on their critical behavior, PRL 94, 218702, 2005.
[Qt] Qt is cross-platform application framework that is widely used for developing
application software with graphical user interface (GUI), qt-project.org
Department of Mathematics, University of Nebraska at Omaha, Omaha, NE 68182
E-mail address: [email protected]
16
APPENDIX
The following is a collection of all the classes used to create this program.
Filename: BinaryTree.hpp
#ifndef BINARYTREE_HPP
#define BINARYTREE_HPP
#include "BinaryTreeNode.h"
#include "defines.h"
template <class T>
class QVector;
class BinaryTree
{
public:
/* Constructor
*/
BinaryTree();
/* Destructor
*/
~BinaryTree();
/* Inserts an item into the binary tree.
* @param[in] aValue item to insert.
* @return success or failure status of the function call.
*/
BNReturn insert(double aValue);
/* Safely destroys the tree.
*/
void destroyTree();
/* Converts the binary tree into a QVector.
* @param[in] aVector vector to place binary tree into.
*/
void treeToVine(QVector<double> &aVector);
int length;
private:
/* Partner to the public method insert.
* @param[in] aValue item to insert.
* @param[in] apNode parent node.
* @return success or failure status of the function call.
*/
BNReturn insert(double aValue, TreeNode* apNode);
/* Partner to the public method destroyTree.
* @param[in] apNode parent node.
*/
void destroyTree(TreeNode* apNode);
/* Partner to the public method treeToVine.
* @param[in] aVector vector to place binary tree into.
* @param[in] apNode node to insert into array.
*/
void treeToVine(QVector<double> &aVector, TreeNode* apNode);
/* Root of the binary tree. */
TreeNode* root;
/* Counter used to convert binary tree to indexed array. */
int counter;
};
#endif // BINARYTREE_HPP
Filename: BinaryTree.cpp
Class Purpose: A binary tree is used a.) when randomly selecting inputs for a node because inserting into the binary tree
structure is the fastest way to find a duplicate entry b.) when searching for cycles (for the same reason).
#include <QVector>
#include "BinaryTree.h"
BinaryTree::BinaryTree() : length(0)
{
root = NULL;
}
BinaryTree::~BinaryTree()
{
destroyTree();
}
void BinaryTree::destroyTree()
{
17
return destroyTree(root);
}
void BinaryTree::destroyTree(TreeNode* apNode)
{
if (apNode != NULL)
{
destroyTree(apNode->left);
destroyTree(apNode->right);
delete apNode;
}
}
BNReturn BinaryTree::insert(double aValue)
{
if (root != NULL)
return insert(aValue, root);
else
{
root = new TreeNode;
root->value = aValue;
root->left = NULL;
root->right = NULL;
length++;
return BN_SUCCESS;
}
}
BNReturn BinaryTree::insert(double aValue, TreeNode* apNode)
{
if (aValue < apNode->value)
{
if (apNode->left != NULL)
return insert(aValue, apNode->left);
else
{
apNode->left = new TreeNode;
apNode->left->value = aValue;
apNode->left->left = NULL;
apNode->left->right = NULL;
length++;
return BN_SUCCESS;
}
}
else if (aValue > apNode->value)
{
if (apNode->right != NULL)
return insert(aValue, apNode->right);
else
{
apNode->right = new TreeNode;
apNode->right->value = aValue;
apNode->right->left = NULL;
apNode->right->right = NULL;
length++;
return BN_SUCCESS;
}
}
else
return BN_ERROR_DUPLICATE;
}
void BinaryTree::treeToVine(QVector<double> &aVector)
{
if (root != NULL)
{
counter = 0;
treeToVine(aVector, root);
}
}
void BinaryTree::treeToVine(QVector<double> &aVector, TreeNode* apNode)
{
aVector[counter] = apNode->value;
counter++;
if (apNode->left != NULL)
treeToVine(aVector, apNode->left);
if (apNode->right != NULL)
treeToVine(aVector, apNode->right);
}
Filename: BinaryTreeNode.hpp
#ifndef BINARYTREENODE_HPP
#define BINARYTREENODE_HPP
18
/*
* Node to a binary tree.
*/
struct TreeNode
{
double value;
TreeNode* left;
TreeNode* right;
};
#endif // BINARYTREENODE_HPP
Filename: TabDialog.hpp
Class Purpose: This is a class that owns the GUI tabs ConfigTab and SimulateTab.
#ifndef TABDIALOG_HPP
#define TABDIALOG_HPP
#include <QDialog>
#include "BooleanNetwork.h"
class QTabWidget;
class ConfigTab;
class SimulateTab;
class TabDialog : public QDialog
{
Q_OBJECT
public:
/* Constructor.
*/
TabDialog(QWidget* parent = 0, BooleanNetwork* bn = 0);
public slots:
/* Shows the simulate tab when an acceptable configuration
* has been inputted and the "Generate Network" button has
* been pressed.
*/
void showSimulateTab();
private:
QTabWidget* mpTabWidget;
ConfigTab* mpConfigWidget;
SimulateTab* mpSimulateWidget;
};
#endif // TABDIALOG_HPP
Filename: TabDialog.cpp
#include
#include
#include
#include
#include
#include
#include
"BooleanNetwork.h"
"ConfigTab.h"
"SimulateTab.h"
"TabDialog.h"
<QtGui>
<QTabWidget>
<QVBoxLayout>
TabDialog::TabDialog(QWidget* parent, BooleanNetwork* bn) :
QDialog(parent),
mpTabWidget(new QTabWidget),
mpConfigWidget(new ConfigTab(0, bn)),
mpSimulateWidget(new SimulateTab(0, bn))
{
mpTabWidget->addTab(mpConfigWidget, tr("Configure"));
QVBoxLayout* mainLayout = new QVBoxLayout;
mainLayout->setSizeConstraint(QLayout::SetNoConstraint);
mainLayout->addWidget(mpTabWidget);
setLayout(mainLayout);
setWindowTitle(tr("MathBio Boolean Network Simulator"));
/* once generateNetwork emits, simulateTab updates its
* parameters */
connect(mpConfigWidget, SIGNAL(networkGenerated()),
mpSimulateWidget, SLOT(updateParameters()));
/* once generateNetwork emits, simulateTab becomes the
* visible tab */
connect(mpSimulateWidget, SIGNAL(readyForDisplay()), this,
SLOT(showSimulateTab()));
}
void TabDialog::showSimulateTab()
{
mpTabWidget->insertTab(1, mpSimulateWidget, tr("Simulate"));
mpTabWidget->setCurrentIndex(1);
}
Filename: ConfigTab.hpp
19
Class Purpose: This class is the GUI tab Configure.
#ifndef CONFIGTAB_HPP
#define CONFIGTAB_HPP
#include <QWidget>
class BooleanNetwork;
namespace Ui
{
class ConfigTab;
}
class ConfigTab : public QWidget
{
Q_OBJECT
public:
/* Constructor.
*/
explicit ConfigTab(QWidget* parent = 0, BooleanNetwork* b = 0);
/* Destructor.
*/
~ConfigTab();
signals:
/* Emits when the boolean network has been successfully generated.
*/
void networkGenerated();
private slots:
/* Slot enables/disables the ability to enter a number into the
* numNodesSpinbox depending on numNodesOption.
*/
void onNumNodesOptionChanged();
/* Slot enables/disables the ability to enter a number into the
* numInputNodesSpinbox depending on numInputNodesOption.
*/
void onNumInputNodesOptionChanged();
/* Slot updates the total percentage label if any of the inputting
* percentages are altered.
*/
void updateTotalPercentage();
/* Either rejects configuration or generates the boolean network.
*/
void on_generatePushbutton_released();
private:
Ui::ConfigTab* ui;
QScopedPointer<BooleanNetwork> mpBn;
};
#endif // CONFIGTAB_HPP
Filename: ConfigTab.cpp
#include
#include
#include
#include
<QMessageBox>
"BooleanNetwork.h"
"ConfigTab.h"
"ui_ConfigTab.h"
ConfigTab::ConfigTab(QWidget *parent, BooleanNetwork* b) :
QWidget(parent),
ui(new Ui::ConfigTab),
mpBn(b)
{
ui->setupUi(this);
bool ok;
ok = connect(ui->randomNumNodes, SIGNAL(clicked()), this,
SLOT(onNumNodesOptionChanged()));
Q_ASSERT(ok);
ok = connect(ui->manualNumNodes, SIGNAL(clicked()), this,
SLOT(onNumNodesOptionChanged()));
Q_ASSERT(ok);
ok = connect(ui->randomNumNodeInputs, SIGNAL(clicked()), this,
SLOT(onNumInputNodesOptionChanged()));
Q_ASSERT(ok);
ok = connect(ui->manualNumNodeInputs, SIGNAL(clicked()), this,
SLOT(onNumInputNodesOptionChanged()));
Q_ASSERT(ok);
ok = connect(ui->thresholdPercentage, SIGNAL(valueChanged(int)),
this, SLOT(updateTotalPercentage()));
Q_ASSERT(ok);
20
ok = connect(ui->canalizing1Percentage,
SIGNAL(valueChanged(int)), this,
SLOT(updateTotalPercentage()));
Q_ASSERT(ok);
ok = connect(ui->canalizing2Percentage,
SIGNAL(valueChanged(int)), this,
SLOT(updateTotalPercentage()));
Q_ASSERT(ok);
ok = connect(ui->biasPercentage, SIGNAL(valueChanged(int)), this,
SLOT(updateTotalPercentage()));
Q_ASSERT(ok); }
ConfigTab::~ConfigTab()
{
delete ui;
}
void ConfigTab::onNumNodesOptionChanged()
{
/* disable numNodesSpinbox if randomNumNodes is selected but
* enable numNodesSpinbox if manualNumNodes is selected */
if (ui->randomNumNodes->isChecked())
ui->numNodesSpinbox->setEnabled(false);
else
ui->numNodesSpinbox->setEnabled(true);
}
void ConfigTab::onNumInputNodesOptionChanged()
{
/* disable numInputNodesSpinbox if randomNumInputs is selected but
* enable numInputNodesSpinbox if manualNumInputs is selected */
if (ui->randomNumNodeInputs->isChecked())
ui->numNodeInputsSpinbox->setEnabled(false);
else
ui->numNodeInputsSpinbox->setEnabled(true);
}
void ConfigTab::updateTotalPercentage()
{
int total = ui->thresholdPercentage->value() +
ui->canalizing1Percentage->value() +
ui->canalizing2Percentage->value() +
ui->biasPercentage->value();
ui->totalPercentage->setValue(total);
}
void ConfigTab::on_generatePushbutton_released()
{
/* generate network if total percentage is 100 percent */
if (ui->totalPercentage->value() != 100)
{
QMessageBox::warning(0, "Error", "Network cannot be generated
unless the total function type percentage is exactly 100
percent", QMessageBox::Ok);
}
if (ui->lowerThreshold->value() >= ui->upperThreshold->value())
{
QMessageBox::warning(0, "Error", "Lower threshold must be less
than upper threshold.", QMessageBox::Ok);
}
if (ui->totalPercentage->value() == 100 &&
ui->lowerThreshold->value() < ui->upperThreshold->value())
{
/* set booleanNetwork's numNodesOption to random or static */
if (ui->randomNumNodes->isChecked())
mpBn->mNumNodesOption = BN_RANDOM;
else
{
mpBn->mNumNodesOption = BN_STATIC;
mpBn->mNumNodes = ui->numNodesSpinbox->value();
}
/* set booleanNetwork's numNodeInputsOption to random or
static */
if (ui->randomNumNodeInputs->isChecked())
mpBn->mNumInputsOption = BN_RANDOM;
else
{
mpBn->mNumInputsOption = BN_STATIC;
mpBn->mNumInputs = ui->numNodeInputsSpinbox->value();
}
/* set booleanNetwork's connectionOption to average or
absolute */
if (ui->average->isChecked())
mpBn->mConnOption = BN_AVERAGE;
else
mpBn->mConnOption = BN_EXACT;
21
/* set threshold function's parameters */
mpBn->mThresholdPercentage = ui->thresholdPercentage->value();
mpBn->mLowerThreshold = ui->lowerThreshold->value() * 100;
mpBn->mUpperThreshold = ui->upperThreshold->value() * 100;
/* set canalizing1 function's parameters */
mpBn->mCanalizing1Percentage = ui->canalizing1Percentage
->value();
mpBn->mRho1 = ui->rho1->value() * 100;
mpBn->mRho2 = ui->rho2->value() * 100;
mpBn->mEta = ui->eta->value() * 100;
/* set canalizing2 function's parameters */
mpBn->mCanalizing2Percentage = ui->canalizing2Percentage
->value();
mpBn->mCanalizing2Bias = ui->canalizing2Bias->value() * 100;
/* set bias function's parameters */
mpBn->mBiasPercentage = ui->biasPercentage->value();
mpBn->mBiasThreshold = ui->biasBias->value() * 100;
if (mpBn->generateNetwork() == BN_SUCCESS)
emit networkGenerated();
}
}
Filename: SimulateTab.hpp
Class Purpose: This class is the GUI tab Simulate
#ifndef SIMULATETAB_HPP
#define SIMULATETAB_HPP
#include <QWidget>
class BooleanNetwork;
namespace Ui
{
class SimulateTab;
}
class SimulateTab : public QWidget
{
Q_OBJECT
public:
/* Constructor.
*/
explicit SimulateTab(QWidget* parent = 0, BooleanNetwork* b = 0);
/* Destructor.
*/
~SimulateTab();
signals:
/* Displays the simulate tab once the network's parameters have
* been updated.
*/
void readyForDisplay();
public slots:
/* Updates the read-only values of the boolean network on the
* simulate page then switches from the configTab to the simulateTab.
*/
void updateParameters();
private slots:
/* Changes T, the amount of cycles the network is simulated over.
*/
void on_periodsSpinbox_valueChanged(int arg1);
/* Starts the network simulation.
*/
void on_simulatePushButton_released();
private:
Ui::SimulateTab* ui;
QScopedPointer<BooleanNetwork> mpBn;
};
#endif // SIMULATETAB_HPP
Filename: SimulateTab.cpp
#include "BooleanNetwork.h"
#include "SimulateTab.h"
#include "ui_SimulateTab.h"
SimulateTab::SimulateTab(QWidget* parent, BooleanNetwork* b) :
QWidget(parent),
ui(new Ui::SimulateTab),
mpBn(b)
{
ui->setupUi(this);
22
}
SimulateTab::~SimulateTab()
{
delete ui;
}
void SimulateTab::updateParameters()
{
/* set the information on the simulate tab to the network's
parameters */
ui->numNodesValue->setText(QString("%0").arg(mpBn->mNumNodes));
ui->numInputsValue->setText(QString("%0").arg(mpBn->mNumInputs));
ui->thresholdValue->setText(QString("%0").arg(mpBn
->mThresholdPercentage));
ui->canalyzing1Value->setText(QString("%0").arg(mpBn
->mCanalizing1Percentage));
ui->canalyzing2Value->setText(QString("%0").arg(mpBn
->mCanalizing2Percentage));
ui->biasValue->setText(QString("%0").arg(mpBn->mBiasPercentage));
ui->averageAbsoluteValue->setText(mpBn->mConnOption == 0?
QString("(Average)") : QString("(Exactly)"));
emit readyForDisplay();
}
void SimulateTab::on_periodsSpinbox_valueChanged(int arg1)
{
mpBn->T = arg1;
}
void SimulateTab::on_simulatePushButton_released()
{
ui->textBrowser->append("Simulation Started.");
ui->textBrowser->append("Working...");
ui->textBrowser->repaint();
BNReturn ret = mpBn->simulateNetwork();
if (ret == BN_SUCCESS)
ui->textBrowser->append("Simulation Complete.");
else
ui->textBrowser->append("Error in Simulation.\n");
}
Filename: defines.hpp
Purpose: This header holds definitions for BooleanNetwork class.
#ifndef DEFINES_HPP
#define DEFINES_HPP
#include <cstdlib>
#include <exception>
#define NOT_NULL(x) if (!x) {return BN_ERROR_NULL_POINTER;}
#define CATCH_ALL catch (std::exception&) {return BN_ERROR;}
/*
* Describes whether a number is to be x or randomly generated.
*/
typedef enum
{
BN_STATIC,
BN_RANDOM
} BNNumberOption;
/*
* Describes whether all nodes in the network have exactly x
* or an average of x inputs.
*/
typedef enum
{
BN_AVERAGE,
BN_EXACT
} BNConnectionOption;
/*
* Possible function types for each node.
*/
typedef enum
{
BN_THRESHOLD,
BN_CANALIZING1,
BN_CANALIZING2,
BN_BIAS
} BNFunctionType;
/*
* Types of canalizing 1 functions.
*/
typedef enum
{
T_POS,
T_NEG,
F_POS,
F_NEG
} BNCan1FunctionTypes;
23
/*
* Types of canalizing 2 functions.
*/
typedef enum
{
T_T_POS,
T_T_NEG,
F_T_POS,
F_T_NEG,
T_F_POS,
T_F_NEG,
F_F_POS,
F_F_NEG
} BNCan2FunctionTypes;
/*
* Possible return values for a function.
*/
typedef enum
{
BN_SUCCESS,
BN_ERROR_NULL_POINTER,
BN_ERROR_PARAM_OUT_OF_RANGE,
BN_ERROR_DUPLICATE,
BN_ERROR,
BN_CYCLE_FOUND,
BN_CYCLE_NOT_FOUND
} BNReturn;
/*
* Possible pairs of parameters to vary.
*/
typedef enum
{
THRESHOLD_LOWER_THRESHOLD_UPPER,
RHO1_RHO2,
PARAM_PAIR_COUNT
} BNParamPairs;
#endif // DEFINES_HPP
Filename: BooleanNetworkNode.hpp
Purpose: Defines a node in the network.
#ifndef BOOLEANNETWORKNODE_H
#define BOOLEANNETWORKNODE_H
#include "BinaryTree.h"
#include "defines.h"
/*
* Structure describing a node in the network.
*/
typedef struct
{
/* State of the node (0 inactive, 1 active). */
int state;
/* Number of input nodes connected to this node. */
int numInputs;
/* Array that stores this node's input nodes indeces in the network.
* An array is useful when randomly selecting indeces.*/
QVector<double> inputs;
/* Node's function type. */
BNFunctionType functionType;
/* If node has type canalizing1, there are four types of canalizing 1
functions. */
BNCan1FunctionTypes can1FunctionType;
/* If node has type canalizing2, there are eight types of canalizing 2
functions. */
BNCan2FunctionTypes can2FunctionType;
/* Index of input node that is canalizing. */
int can1Index;
/* Index of input node that is canalizing. */
int can2Index;
} BNNode;
#endif // BOOLEANNETWORKNODE_HPP
Filename: BooleanNetwork.hpp
Class Purpose: This class generates and simulates the Boolean network.
#ifndef BOOLEANNETWORK_H
#define BOOLEANNETWORK_H
#include <fstream>
24
#include
#include
#include
#include
<iostream>
<QVector>
"BooleanNetworkNode.h"
"defines.h"
class BooleanNetwork
{
public:
/*
* Constructor.
*/
BooleanNetwork();
/*
* Destructor.
*/
~BooleanNetwork();
/* If the check box for a param pair in the GUI gets checked, this
* function is called.
* @return success or failure status for the function call.
*/
/* Generate the boolean network that will be used for simulation.
* @return success or failure status of the function call.
*/
BNReturn generateNetwork();
/* Runs the network through to time T then checks the stability of the
* network by examining each node.
* @return success or failure status of the function call.
*/
BNReturn simulateNetwork();
int mNumNodes;
int mLowerNumNodes;
int mUpperNumNodes;
BNNumberOption mNumNodesOption;
int mNumInputs;
int mLowerNumInputs;
int mUpperNumInputs;
BNNumberOption mNumInputsOption;
BNConnectionOption mConnOption;
int
int
int
int
int
mThresholdPercentage;
mCanalizing1Percentage;
mCanalizing2Percentage;
mBiasPercentage;
mTotalPercentage;
int mRho1;
int mRho2;
int mEta;
/* Fraction of canalyzing 1 functions
with canalized output */
/* Bias of function OTHER. */
/* Fraction of canalyzing 1 functions
with canalized intput */
int mCanalizing2Bias;
int mLowerThreshold;
int mUpperThreshold;
int mBiasThreshold;
int T;
private:
/* Flips each node in the network sequentially and determines if the
* network recovered.
* @return stability of the network.
*/
double flipNodes();
/* Evaluates each node at time t until time T is reached or a cycle is
* found.
* @param tree binary tree that stores cycle.
* @return success or failure status of the function call.
*/
BNReturn evaluateToTimeT(BinaryTree &tree);
/* Evaluates a node at time t with a threshold function.
* @param[in] aNodeIndex index of node being evaluated.
* @param[out] aNewState new state of booleanNetworkNode (0 or 1).
* @return success or failure status of the function call.
*/
BNReturn evalThresholdNode(int aNodeIndex, int &aNewState);
/* Evaluates a node at time t with a canalyzing1 function.
* @param[in] aNodeIndex index of node being evaluated.
* @param[out] aNewState new state of booleanNetworkNode (0 or 1).
* @return success or failure status of the function call.
*/
BNReturn evalCanalizing1Node(int aNodeIndex, int &aNewState);
/* Evaluates a node at time t with a canalyzing2 function.
* @param[in] aNodeIndex index of node being evaluated.
25
* @param[out] aNewState new state of booleanNetworkNode (0 or 1).
* @return success or failure status of the function call.
*/
BNReturn evalCanalizing2Node(int aNodeIndex, int &aNewState);
/* Evaluates a node at time t with a bias function.
* @param[out] aNewState new state of booleanNetworkNode (0 or 1).
* @return success or failure status of the function call.
*/
BNReturn evalBiasNode(int &aNewState);
/* Checks network state to see if a cycle has been reached.
* @param[in] tree binary tree storing network's state in the form of
* a decimal.
* @return success or failure status of the function call.
*/
BNReturn checkForCycle(BinaryTree &tree);
/* Varys parameter pairs.
* @param[in] x value for param 1.
* @param[in] y value for param 2.
*/
void varyParamPair(int pair, double x, double y);
/* Returns the appropriate filename for a varied parameter pair.
* @param[in] param pair to get file name for.
* @return file name.
*/
QString getFileName(int paramPair);
/* Assigns each canalyzing1 node a canalyzing1 node function in
* accordance with Rho1 and Eta.
* @return success or failure status of the function call.
*/
BNReturn assignCanalizing1Functions();
/* Assigns each canalizing2 node a canalizing2 node function.
* @return success or failure status of the function call.
*/
BNReturn assignCanalizing2Functions();
QVector<BNNode> mNetwork;
QVector<bool> mParamPairs;
/* for testing purposes only */
std::ofstream networkInfo;
};
#endif // BOOLEANNETWORK_HPP
Filename: BooleanNetwork.cpp
#include
#include
#include
#include
#include
<cmath>
<cstdlib>
<time.h>
<QString>
"BooleanNetwork.h"
//#define TESTING_ENABLED
#define NETWORK_RECOVERED
#define NETWORK_NOT_RECOVERED
/* for rand() and srand() */
/* for time() and typedef time_t */
1
0
BooleanNetwork::BooleanNetwork() :
mNumNodes(5),
mLowerNumNodes(3),
/* number of nodes in network */
/* lower bound for mNumNodes if
mNumNodes is random */
mUpperNumNodes(25),
/* upper bound for mNumNodes if
mNumNodesOption(BN_STATIC),
/* determines whether mNumNodes =
mNumNodes is random */
mNumNodes or if mNumNodes is
randomly generated in range
lowermNumNodes to uppermNumNodes
*/
mNumInputs(2),
mLowerNumInputs(1),
/* number of input nodes per node */
/* lower bound for numInputs if
mUpperNumInputs(3),
/* upper bound for numInputs if
mNumInputsOption(BN_STATIC),
/* determines whether numInputs =
ConnOption = BN_AVERAGE */
ConnOption = BN_AVERAGE */
mConnOption(BN_EXACT),
numInputs or if numInputs is
randomly generated in range
lowerNumInputs to upperNumInputs
*/
/* determines whether each node has
exactly numInputs or an
average of numInputs across
network */
mThresholdPercentage(25),
/* percentage of nodes w/threshold
mCanalizing1Percentage(25),
/* percentage of nodes w/canalyzing
function */
function w/1 input */
26
mCanalizing2Percentage(25),
/* percentage of nodes w/canalyzing
function w/2 inputs */
/* percentage of nodes w/bias function
*/
/* sum of previous four parameters */
mBiasPercentage(25),
mTotalPercentage(100),
mRho1(50),
/* percentage of canalyzing1 functions
with canalized output */
/* bias of OTHER function */
/* percentage of canalyzing1 functions
with canalized input */
mRho2(50),
mEta(10),
mLowerThreshold(30),
/* lower boundary for threshold
mUpperThreshold(60),
/* upper boundary for threshold
mBiasThreshold(50),
/* probability of node being on for
T(100),
/* cycles to simulate over */
function */
function */
bias function */
mNetwork(),
mParamPairs(PARAM_PAIR_COUNT, false)
{
#ifdef TESTING_ENABLED
networkInfo.open("networkInfoForTestingPurposes.txt");
networkInfo << "Boolean Network Simulation (For Testing Purposes)\n\n";
#endif
}
BooleanNetwork::~BooleanNetwork()
{
networkInfo.close();
}
BNReturn BooleanNetwork::generateNetwork()
{
try
{
/* prepare for random number generation */
time_t seconds;
time(&seconds);
srand((unsigned int) seconds);
int lower = 1;
int upper = 1;
/* check to see if the number of nodes for the network was
* initialized by the user or if it needs to be randomly generated
*/
if (mNumNodesOption == BN_RANDOM)
mNumNodes = rand() % (mUpperNumNodes - mLowerNumNodes + 1) +
mLowerNumNodes;
/* create the network */
mNetwork = QVector<BNNode>(mNumNodes);
#ifdef TESTING_ENABLED
networkInfo << "SETUP\n";
if (mNumNodesOption)
{
networkInfo << "Num Nodes
networkInfo << "Upper Num
networkInfo << "Lower Num
}
else
networkInfo << "Num Nodes
Option:\tBN_RANDOM\n";
Nodes:\t" << mUpperNumNodes << "\n";
Nodes:\t" << mLowerNumNodes << "\n";
Option:\tBN_STATIC\n";
networkInfo << "Num Nodes:\t\t" << mNumNodes << "\n";
if (mConnOption)
networkInfo << "Num
"
else
networkInfo << "Num
"
Inputs:\t\t" << mNumInputs <<
(absolute)\n";
Inputs:\t\t" << mNumInputs <<
(average)\n";
networkInfo << "Node Function Type\n";
networkInfo << "\tThreshold:\t" << mThresholdPercentage << "\n";
networkInfo << "\tCanalyzing 1:\t" << mCanalizing1Percentage <<
"\n";
networkInfo << "\tCanalyzing 2:\t" << mCanalizing2Percentage <<
"\n";
networkInfo << "\tBias:\t\t" << mBiasPercentage << "\n\n";
#endif
/* prepare for function type assignment */
int thresholdLower = 0;
int thresholdUpper = 0;
int canalizing1Lower = 0;
int canalizing1Upper = 0;
int canalizing2Lower = 0;
int canalizing2Upper = 0;
int biasLower = 0;
int biasUpper = 0;
27
thresholdUpper = mThresholdPercentage;
canalizing1Lower = thresholdUpper;
if (mThresholdPercentage)
canalizing1Lower++;
canalizing1Upper = thresholdUpper + mCanalizing1Percentage;
canalizing2Lower = canalizing1Upper;
if (mCanalizing1Percentage)
canalizing2Lower++;
canalizing2Upper = canalizing1Upper + mCanalizing2Percentage;
biasLower = canalizing2Upper;
if (mCanalizing2Percentage)
biasLower++;
biasUpper = canalizing2Upper + mBiasPercentage;
/* initialize the parameters for each node of the network */
for (int n=0; n<mNumNodes; n++)
{
/* STATE : randomly assigned during simulation */
/* NUMBER OF INPUTS */
int numInputs = -1;
if (mConnOption == BN_EXACT)
numInputs = mNumInputs;
else
{
lower = 1;
upper = mNumInputs + (mNumInputs - lower);
numInputs = rand() % (upper - lower + 1) + lower;
}
mNetwork[n].numInputs = numInputs;
mNetwork[n].inputs = QVector<double>(numInputs);
/* INPUT NODES : create a binary tree that stores each input
* node's network index */
lower = 0;
upper = mNumNodes-1;
BinaryTree tree;
for (int i=0; i<mNetwork[n].numInputs; i++)
{
/* try picking a random input node. reject it if the node
* inputs to itself or if that node has already been
* selected as an input node */
BNReturn ret = BN_ERROR;
while (ret != BN_SUCCESS)
{
/* pick random node */
int nodeIndex = rand() % (upper - lower + 1) + lower;
/* check to see if input node is the node itself */
if (nodeIndex == n)
ret = BN_ERROR;
/* try inputing the input node into this node's binary
* tree of input nodes. if it is not successfully
* inputed, that means the node
* has already been selected as an input */
ret = tree.insert(nodeIndex);
}
}
/* initialize mNetwork[n].inputs to values stored in binary
tree */
tree.treeToVine(mNetwork[n].inputs);
/* FUNCTION TYPE */
lower = 1;
upper = 100;
int r = rand() % (upper - lower + 1) + lower;
if (mNetwork[n].numInputs == 1)
{
while (r >= canalizing2Lower && r <= canalizing2Upper)
{
r = rand() % (upper - lower + 1) + lower;
}
}
BNFunctionType f;
if (r >= thresholdLower && r <= thresholdUpper)
f = BN_THRESHOLD;
else if (r >= canalizing1Lower && r <= canalizing2Upper)
f = BN_CANALIZING1;
else if (r >= canalizing2Lower && r <= canalizing2Upper)
f = BN_CANALIZING2;
else if (r >= biasLower && r <= biasUpper)
f = BN_BIAS;
mNetwork[n].functionType = f;
/* CANALIZING VALUES */
mNetwork[n].can1Index = -1;
mNetwork[n].can2Index = -1;
28
}
/* if a function is of type canalyzing1, it needs to be assigned
* specific canalyzing1 function */
this->assignCanalizing1Functions();
}
CATCH_ALL
return BN_SUCCESS;
}
BNReturn BooleanNetwork::simulateNetwork()
{
try
{
mParamPairs[RHO1_RHO2] = true;
for (int paramPair=0; paramPair<PARAM_PAIR_COUNT; paramPair++)
{
if (mParamPairs[paramPair])
{
std::ofstream file;
/* ridiculous work around to get QString filename into
correct format */
QString a = this->getFileName(paramPair);
QByteArray b = a.toLatin1();
const char* filename = b.data();
file.open(filename);
file
file
file
file
file
<<
<<
<<
<<
<<
file <<
file <<
file <<
file <<
file <<
"%X Intervals: 30\n";
"%Y Intervals: 30\n";
"%Number of Nodes: " << mNumNodes << "\n";
"%Connectivity: " << mNumInputs << "\n";
"%Percent Threshold: " << mThresholdPercentage <<
"\n";
"%Percent Canalyzing 1: " <<
mCanalizing1Percentage << "\n";
"%Percent Canalyzing 2: " <<
mCanalizing2Percentage << "\n";
"%Percent Bias: " << mBiasPercentage << "\n";
"%T: " << T << "\n";
"% Eta: " << mEta / 100.0 << "\n\n";
//this->flipNodes();
double stability;
for (double x=0; x<=30; x++)
{
for (double y=0; y<=30; y++)
{
std::cout << "x = " << x << ", y = " << y <<
std::endl;
this->varyParamPair(paramPair, x, y);
stability = this->flipNodes();
file << stability << "\t";
std::cout << "stability = " << stability <<
std::endl;
}
file << "\n";
}
file.close();
}
}
/* create coordinates files */
std::ofstream coordinatesX;
std::ofstream coordinatesY;
coordinatesX.open("coordinatesX.txt");
coordinatesY.open("coordinatesY.txt");
for (double lower=0; lower<30; lower++)
{
for (double upper=0; upper<=30; upper++)
{
coordinatesX << (upper/30) << "\t";
coordinatesY << (lower/30) << "\t";
}
coordinatesX << "\n";
coordinatesY << "\n";
}
coordinatesX.close();
coordinatesY.close();
}
CATCH_ALL
return BN_SUCCESS;
}
double BooleanNetwork::flipNodes()
{
QVector<int> impact(mNumNodes, 0);
for (int node=0; node<mNumNodes; node++)
29
{
for (int times=0; times<50; times++)
{
/* randomly initialize each node with 0 or 1 */
for (int i=0; i<mNumNodes; i++)
mNetwork[i].state = rand() % 2;
/* binary tree stores the network state by condensing the
* sequence of 0s and 1s into a decimal number */
BinaryTree tree;
/* simulate network until a cycle is reached (either an actual
* cycle or the cycle of t=0 to t=T) */
this->evaluateToTimeT(tree);
/* flip node */
mNetwork[node].state = !(mNetwork[node].state);
/* binary tree stores the network state by condensing the
* sequence of 0s and 1s into a decimal number */
BinaryTree tree2;
/* simulate network and search for a cycle. did network
recover? */
this->evaluateToTimeT(tree2);
/* compare trees to see if the network reached the same state
as before */
QVector<double> network(tree.length);
QVector<double> network2(tree2.length);
tree.treeToVine(network);
tree2.treeToVine(network2);
bool match = false;
int size = network.size();
int size2 = network2.size();
for (int i=0; i<size; i++)
{
for (int j=0; j<size2; j++)
{
if (network[i] == network2[j])
{
match = true;
break;
}
}
if (match == true)
break;
}
if (match)
impact[node] += NETWORK_RECOVERED;
else
impact[node] += NETWORK_NOT_RECOVERED;
}
}
/* average the stability over each node */
double stability = 0;
for (int node=0; node<mNumNodes; node++)
{
stability += impact[node] / 50.0;
}
/* use the average stability of each node to create a network average
*/
return stability = stability / mNumNodes;
}
BNReturn BooleanNetwork::evaluateToTimeT(BinaryTree &tree)
{
BNReturn ret = BN_ERROR;
try
{
/* simulate the network until time T or a cycle has been reached
*/
for (int t=0; t<T; t++)
{
#ifdef TESTING_ENABLED
networkInfo << "SIMULATION (t=" << t << ")\n\n";
networkInfo << "State: ";
for (int node=0; node<mNumNodes; node++)
networkInfo << mNetwork[node].state << " ";
networkInfo << "\n";
#endif
/* check to see if a cycle has been reached...if it has, we're
done */
ret = this->checkForCycle(tree);
if (ret == BN_CYCLE_FOUND)
{
#ifdef TESTING_ENABLED
30
networkInfo << "Cycle Found !!!\n\n";
#endif
break;
}
else
{
/* create an QVector to store the state of each node at
time t */
QVector<int> newNetworkState(mNumNodes);
/* evaluate the new state of each node in the network */
for (int node=0; node<mNumNodes; node++)
{
#ifdef TESTING_ENABLED
networkInfo << "Node: " << node << "\n";
networkInfo << "State: " << mNetwork[node].state <<
"\n";
networkInfo << "Type: ";
if (mNetwork[node].functionType == BN_THRESHOLD)
networkInfo << "Threshold\n";
else if (mNetwork[node].functionType ==
BN_CANALIZING1)
networkInfo << "Canalyzing 1\n";
else if (mNetwork[node].functionType ==
BN_CANALIZING2)
networkInfo << "Canalyzing 2\n";
else if (mNetwork[node].functionType == BN_BIAS)
networkInfo << "Bias\n";
#endif
/* evaluate state */
int newState = -1;
if (mNetwork[node].functionType == BN_THRESHOLD)
{
ret = evalThresholdNode(node, newState);
Q_ASSERT(ret == BN_SUCCESS);
}
else if (mNetwork[node].functionType ==
BN_CANALIZING1)
{
ret = evalCanalizing1Node(node, newState);
Q_ASSERT(ret == BN_SUCCESS);
}
else if (mNetwork[node].functionType ==
BN_CANALIZING2)
{
ret = evalCanalizing2Node(node, newState);
Q_ASSERT(ret == BN_SUCCESS);
}
else
{
Q_ASSERT(mNetwork[node].functionType == BN_BIAS);
ret = evalBiasNode(newState);
Q_ASSERT(ret == BN_SUCCESS);
}
newNetworkState[node] = newState;
}
for (int node=0; node<mNumNodes; node++)
mNetwork[node].state = newNetworkState[node];
} /* end else (cycle not found) */
} /* end for (int t=0; t<T; t++) */
}
CATCH_ALL
return ret;
}
BNReturn BooleanNetwork::evalThresholdNode(int aNodeIndex, int &aNewState)
{
NOT_NULL(&aNewState);
try
{
double numActive = 0;
double percentActive = 0.0;
/* find percentage of active input nodes */
int size = (mNetwork[aNodeIndex].inputs).size();
for (int node=0; node<size; node++)
{
int inputNode = mNetwork[aNodeIndex].inputs[node];
if (mNetwork[inputNode].state == 1)
numActive++;
}
percentActive = numActive / (mNetwork[aNodeIndex].inputs).size() *
100;
if (percentActive > mUpperThreshold ||
percentActive < mLowerThreshold)
aNewState = 0;
31
else
aNewState = 1;
#ifdef TESTING_ENABLED
networkInfo << "Input Indexes:\n";
for (int i=0; i<(mNetwork[aNodeIndex].inputs).size(); i++)
networkInfo << mNetwork[aNodeIndex].inputs[i] << " ";
networkInfo << "\n";
networkInfo << "Input States:\n";
for (int i=0; i<(mNetwork[aNodeIndex].inputs).size(); i++)
networkInfo << mNetwork[mNetwork[aNodeIndex].inputs[i]].state
<< " ";
networkInfo << "\n";
networkInfo
networkInfo
networkInfo
networkInfo
<<
<<
<<
<<
"Lower Threshold: " << mLowerThreshold << "\n";
"Upper Threshold: " << mUpperThreshold << "\n";
"Percent Active: " << percentActive << "\n";
"New State: " << aNewState << "\n\n";
#endif
}
CATCH_ALL
return BN_SUCCESS;
}
BNReturn BooleanNetwork::evalCanalizing1Node(int aNodeIndex, int &aNewState)
{
NOT_NULL(&aNewState);
try
{
if (mNetwork[aNodeIndex].can1Index == -1)
{
/* select a random input node to be canalizing */
int lower = 0;
int upper = (mNetwork[aNodeIndex].inputs).size() - 1;
int index = rand() % (upper - lower + 1) + lower;
/* translate from local (array of inputs) to global (overall
network array) */
mNetwork[aNodeIndex].can1Index =
mNetwork[aNodeIndex].inputs[index];
}
BNCan1FunctionTypes type = mNetwork[aNodeIndex].can1FunctionType;
int state = mNetwork[mNetwork[aNodeIndex].can1Index].state;
int r = -1;
if (type == T_POS && state == 1)
aNewState = 1;
else if (type == T_NEG && state == 1)
aNewState = 0;
else if (type == F_POS && state == 0)
aNewState = 1;
else if (type == F_NEG && state == 0)
aNewState = 0;
else
{
/* generate random number to check against bias */
int lower = 0;
int upper = 100;
r = rand() % (upper - lower + 1) + lower;
if (r >= mRho2)
aNewState = 1;
else
aNewState = 0;
}
#ifdef TESTING_ENABLED
networkInfo << "Input Indexes:\n";
for (int i=0; i<(mNetwork[aNodeIndex].inputs).size(); i++)
networkInfo << mNetwork[aNodeIndex].inputs[i] << " ";
networkInfo << "\n";
networkInfo << "Input States:\n";
for (int i=0; i<(mNetwork[aNodeIndex].inputs).size(); i++)
networkInfo << mNetwork[mNetwork[aNodeIndex].inputs[i]].state
<< " ";
networkInfo << "\n";
networkInfo << "Canalyzing Input Index: " <<
mNetwork[aNodeIndex].can1Index << "\n";
networkInfo << "Canalyzing Node State: " << state << "\n";
networkInfo << "Random Number Generated: " << r << "\n";
networkInfo << "Rho2 (canalyzing1 bias): " << mRho2 << "\n";
networkInfo << "Canalyzing1 Function Type: ";
if (mNetwork[aNodeIndex].can1FunctionType == T_POS)
networkInfo << "1->1\n";
else if (mNetwork[aNodeIndex].can1FunctionType == T_NEG)
networkInfo << "1->0\n";
else if (mNetwork[aNodeIndex].can1FunctionType == F_POS)
networkInfo << "0->1\n";
32
else if (mNetwork[aNodeIndex].can1FunctionType == F_NEG)
networkInfo << "0->0\n";
networkInfo << "New State: " << aNewState << "\n\n";
#endif
}
CATCH_ALL
return BN_SUCCESS;
}
BNReturn BooleanNetwork::evalCanalizing2Node(int aNodeIndex, int &aNewState)
{
NOT_NULL(&aNewState);
try
{
if (mNetwork[aNodeIndex].can1Index == -1)
{
/* select 2 random canalyzing input nodes */
int lower = 0;
int upper = (mNetwork[aNodeIndex].inputs).size() - 1;
int index = rand() % (upper - lower + 1) + lower;
int index2 = index;
while (index2 == index)
{
index2 = rand() % (upper - lower + 1) + lower;
}
/* translate from local (array of inputs) to global (overall
network array) */
mNetwork[aNodeIndex].can1Index =
mNetwork[aNodeIndex].inputs[index];
mNetwork[aNodeIndex].can2Index =
mNetwork[aNodeIndex].inputs[index2];
}
BNCan2FunctionTypes type = mNetwork[aNodeIndex].can2FunctionType;
int state = mNetwork[mNetwork[aNodeIndex].can1Index].state;
int state2 = mNetwork[mNetwork[aNodeIndex].can2Index].state;
int r = -1;
if (type == T_T_POS && state == 1
aNewState = 1;
else if (type == T_T_NEG && state
aNewState = 0;
else if (type == T_F_POS && state
aNewState = 1;
else if (type == T_F_NEG && state
aNewState = 0;
else if (type == F_T_POS && state
aNewState = 1;
else if (type == F_T_POS && state
aNewState = 0;
else if (type == F_F_POS && state
aNewState = 1;
else if (type == F_F_POS && state
aNewState = 0;
else
{
int lower = 0;
int upper = 100;
r = rand() % (upper - lower +
&& state2 == 1)
== 1 && state2 == 1)
== 1 && state2 == 0)
== 1 && state2 == 0)
== 0 && state2 == 1)
== 0 && state2 == 1)
== 0 && state2 == 0)
== 0 && state2 == 0)
1) + lower;
if (r >= mCanalizing2Bias)
aNewState = 1;
else
aNewState = 0;
}
#ifdef TESTING_ENABLED
networkInfo << "Input Indexes:\n";
for (int i=0; i<(mNetwork[aNodeIndex].inputs).size(); i++)
networkInfo << mNetwork[aNodeIndex].inputs[i] << " ";
networkInfo << "\n";
networkInfo << "Input States:\n";
for (int i=0; i<(mNetwork[aNodeIndex].inputs).size(); i++)
networkInfo << mNetwork[mNetwork[aNodeIndex].inputs[i]].state
<< " ";
networkInfo << "\n";
networkInfo << "Canalyzing Input 1 Index: " <<
mNetwork[aNodeIndex].can1Index << "\n";
networkInfo << "Canalizing Input 2 Index: " <<
mNetwork[aNodeIndex].can2Index << "\n";
networkInfo << "Canalyzing 1 Node State: " << state << "\n";
networkInfo << "Canalizing 2 Node State: " << state2 << "\n";
networkInfo << "Random Number Generated: " << r << "\n";
networkInfo << "Canalizing 2 Bias: " << mRho2 << "\n";
networkInfo << "Canalizing2 Function Type: ";
if (mNetwork[aNodeIndex].can2FunctionType == T_T_POS)
networkInfo << "11->1\n";
else if (mNetwork[aNodeIndex].can2FunctionType == T_T_NEG)
networkInfo << "11->0\n";
33
else if (mNetwork[aNodeIndex].can2FunctionType == F_T_POS)
networkInfo << "01->1\n";
else if (mNetwork[aNodeIndex].can2FunctionType == F_T_NEG)
networkInfo << "01->0\n";
else if (mNetwork[aNodeIndex].can2FunctionType == T_F_POS)
networkInfo << "10->1\n";
else if (mNetwork[aNodeIndex].can2FunctionType == T_F_NEG)
networkInfo << "10->0\n";
else if (mNetwork[aNodeIndex].can2FunctionType == F_F_POS)
networkInfo << "00->1\n";
else if (mNetwork[aNodeIndex].can2FunctionType == F_F_NEG)
networkInfo << "00->0\n";
networkInfo << "New State: " << aNewState << "\n\n";
#endif
}
CATCH_ALL
return BN_SUCCESS;
}
BNReturn BooleanNetwork::evalBiasNode(int &aNewState)
{
NOT_NULL(&aNewState);
try
{
int lower = 0;
int upper = 100;
int r = rand() % (upper - lower + 1) + lower;
if (r >= mBiasThreshold)
aNewState = 1;
else
aNewState = 0;
#ifdef TESTING_ENABLED
networkInfo << "Bias Threshold: " << mBiasThreshold << "\n";
networkInfo << "Random Number: " << r << "\n";
networkInfo << "New State: " << aNewState << "\n\n";
#endif
}
CATCH_ALL
return BN_SUCCESS;
}
BNReturn BooleanNetwork::assignCanalizing1Functions()
{
try
{
int percentA = 0;
/* percent of canalyzing1 functions with the
function CAN_POS */
int percentB = 0;
/* percent of canalyzing1 functions with the
function NOT_CAN_POS */
int percentC = 0;
/* percent of canalyzing1 functions with the
function NOT_CAN_NEG */
if (mRho1 == 0)
percentC = mEta;
else if (mRho1
{
percentA =
percentB =
percentC =
}
else if (mRho1
{
percentA =
percentB =
percentC =
}
int
int
int
int
int
int
int
int
lowerA
upperA
lowerB
upperB
lowerC
upperC
lowerD
upperD
=
=
=
=
=
=
=
=
<= mEta)
mRho1 / 2;
mRho1 / 2;
mEta = mRho1 / 2;
> mEta)
mEta / 2;
mRho1 - mEta / 2;
mEta / 2;
0;
0;
0;
0;
0;
0;
0;
100;
upperA = percentA;
lowerB = upperA;
if (percentA)
lowerB++;
upperB = upperA + percentB;
lowerC = upperB;
if (percentB)
lowerC++;
34
upperC = upperB + percentC;
lowerD = upperC;
if (percentC)
lowerD++;
int lower = 0;
int upper = 100;
int r;
for (int node=0; node<mNumNodes; node++)
{
if (mNetwork[node].functionType == BN_CANALIZING1)
{
/* generate a random number 0-100 */
r = rand() % (upper - lower + 1) + lower;
if (lowerA <= r && r <= upperA)
mNetwork[node].can1FunctionType = T_POS;
else if (lowerB <= r && r <= upperB)
mNetwork[node].can1FunctionType = F_POS;
else if (lowerC <= r && r <= upperC)
mNetwork[node].can1FunctionType = F_NEG;
else
{
Q_ASSERT(lowerD <= r && r <= upperD);
mNetwork[node].can1FunctionType = T_NEG;
}
}
}
#ifdef TESTING_ENABLED
double rho1 = 0;
double eta = 0;
for (int node=0; node<mNumNodes; node++)
{
if (mNetwork[node].can1FunctionType == T_POS ||
mNetwork[node].can1FunctionType == F_POS)
{
rho1++;
}
if (mNetwork[node].can1FunctionType == T_POS ||
mNetwork[node].can1FunctionType == T_NEG)
{
eta++;
}
}
rho1 = rho1 / mNumNodes;
eta = eta / mNumNodes;
int x = 7;
#endif
}
CATCH_ALL
return BN_SUCCESS;
}
BNReturn BooleanNetwork::assignCanalizing2Functions()
{
try
{
int lower = 1;
int upper = 8;
int r = -1;
for (int node=0; node<mNumNodes; node++)
{
if (mNetwork[node].functionType == BN_CANALIZING2)
{
r = rand() % (upper - lower + 1) + lower;
switch(r)
{
case 1:
mNetwork[node].can2FunctionType = T_T_POS;
break;
case 2:
mNetwork[node].can2FunctionType = T_T_NEG;
break;
case 3:
mNetwork[node].can2FunctionType = T_F_POS;
break;
case 4:
mNetwork[node].can2FunctionType = T_F_NEG;
break;
35
case 5:
mNetwork[node].can2FunctionType = F_T_POS;
break;
case 6:
mNetwork[node].can2FunctionType = F_T_NEG;
break;
case 7:
mNetwork[node].can2FunctionType = F_F_POS;
break;
case 8:
mNetwork[node].can2FunctionType = F_F_NEG;
break;
}
}
}
}
CATCH_ALL
return BN_SUCCESS;
}
BNReturn BooleanNetwork::checkForCycle(BinaryTree &tree)
{
/* convert the network's current state, a list of 0s and 1s, to a
decimal number */
double decimal = 0;
int n = mNumNodes;
for (int node=0; node<mNumNodes; node++)
{
#ifdef TESTING_ENABLED
int state = mNetwork[node].state;
#endif
decimal += mNetwork[node].state * pow(2.0, n-1);
n -= 1;
}
/* insert decimal number into binary tree. if number cannot be
inserted then
* the binary tree already contains that number, meaning a cycle has
been reached */
BNReturn ret = tree.insert(decimal);
return ret == BN_SUCCESS? BN_CYCLE_NOT_FOUND : BN_CYCLE_FOUND;
}
void BooleanNetwork::varyParamPair(int pair, double x, double y)
{
if (pair == THRESHOLD_LOWER_THRESHOLD_UPPER)
{
mLowerThreshold = y/30*100;
mUpperThreshold = x/30*100;
}
else if (pair == RHO1_RHO2)
{
mRho1 = y/30*100;
mRho2 = x/30*100;
}
}
QString BooleanNetwork::getFileName(int paramPair)
{
QString filename = "";
if (paramPair == THRESHOLD_LOWER_THRESHOLD_UPPER)
filename = "(Threshold Lower)(Threshold Upper).txt";
else if (paramPair == RHO1_RHO2)
filename = "(Rho1)(Rho2).txt";
return filename;
}
Filename: main.cpp
Purpose: Runs the program.
#include "TabDialog.h"
#include <QApplication>
int main(int argc, char *argv[])
{
QApplication app(argc, argv);
/* create a boolean network */
BooleanNetwork booleanNetwork;
/* create the MathBio Boolean Network Simulator environment */
TabDialog tabdialog(0, &booleanNetwork);
tabdialog.setMinimumSize(358,458);
tabdialog.show();
36
return app.exec();
}
ORGANIZATION
TabDialog
BooleanNetwork
ConfigTab
SimulateTab
The class TabDialog holds pointers to the classes: ConfigTab, SimulateTab and
BooleanNetwork. The BooleanNetwork pointer is a QScopedPointer so it may be safely
shared with the classes ConfigTab and SimulateTab. ConfigTab uses the pointer to initialize
BooleanNetwork’s parameters and SimulateTab uses pointer to show parameters on
SimulateTab.
Fly UP