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.