class27_booleannetwork_R_template.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class 27 - Boolean Networks\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define the names of the nodes (molecules) in the model\n",
    "\n",
    "nodes <- c(\"Cell Size\",\n",
    "           \"Cln3\",\n",
    "           \"MBF\",\n",
    "           \"Clb5,6\",\n",
    "           \"Mcm1/SFF\",\n",
    "           \"Swi5\",\n",
    "           \"Sic1\",\n",
    "           \"Clb1,2\",\n",
    "           \"Cdc20&Cdc14\",\n",
    "           \"Cdh1\",\n",
    "           \"Cln1,2\",\n",
    "           \"SBF\")\n",
    "           \n",
    "# count the number of nodes in the network\n",
    "N <- length(nodes)\n",
    "\n",
    "# define an NxN matrix of dependencies (these are the \"edges\" in the network)\n",
    "# the row (first index) identifies the \"regulator\", and the column\n",
    "# (second index) identifies the target molecule.\n",
    "# NOTE:  I define aij as the transpose of how Li et al. define theirs!\n",
    "a <- matrix(nrow=N, ncol=N)\n",
    "\n",
    "# label the rows and columns of the matrix by the node names\n",
    "rownames(a) <- nodes\n",
    "colnames(a) <- nodes\n",
    "\n",
    "# initialize all entries of the matrix to zero\n",
    "a[,] <- 0\n",
    "\n",
    "# define the activation edges (green arrows)\n",
    "a[\"Cell Size\",    \"Cln3\"]        <- 1\n",
    "a[\"Cln3\",         \"MBF\"]         <- 1\n",
    "a[\"Cln3\",         \"SBF\"]         <- 1\n",
    "a[\"SBF\",          \"Cln1,2\"]      <- 1\n",
    "a[\"MBF\",          \"Clb5,6\"]      <- 1\n",
    "a[\"Cdc20&Cdc14\",  \"Sic1\"]        <- 1\n",
    "a[\"Clb5,6\",       \"Mcm1/SFF\"]    <- 1\n",
    "a[\"Clb5,6\",       \"Clb1,2\"]      <- 1\n",
    "a[\"Mcm1/SFF\",     \"Clb1,2\"]      <- 1\n",
    "a[\"Clb1,2\",       \"Mcm1/SFF\"]    <- 1\n",
    "a[\"Mcm1/SFF\",     \"Cdc20&Cdc14\"] <- 1\n",
    "a[\"Mcm1/SFF\",     \"Swi5\"]        <- 1\n",
    "a[\"Cdc20&Cdc14\",  \"Swi5\"]        <- 1\n",
    "a[\"Cdc20&Cdc14\",  \"Cdh1\"]        <- 1\n",
    "a[\"Clb1,2\",       \"Cdc20&Cdc14\"] <- 1\n",
    "a[\"Swi5\",         \"Sic1\"]        <- 1\n",
    "\n",
    "# define the repression edges (red arrows)\n",
    "a[\"Clb1,2\",       \"SBF\"]         <- -1\n",
    "a[\"Clb1,2\",       \"MBF\"]         <- -1\n",
    "a[\"Clb1,2\",       \"Cdh1\"]        <- -1\n",
    "a[\"Clb1,2\",       \"Swi5\"]        <- -1\n",
    "a[\"Clb1,2\",       \"Sic1\"]        <- -1\n",
    "a[\"Cln1,2\",       \"Sic1\"]        <- -1\n",
    "a[\"Clb5,6\",       \"Sic1\"]        <- -1\n",
    "a[\"Sic1\",         \"Clb5,6\"]      <- -1\n",
    "a[\"Cln1,2\",       \"Cdh1\"]        <- -1\n",
    "a[\"Clb5,6\",       \"Cdh1\"]        <- -1\n",
    "a[\"Cdc20&Cdc14\",  \"Clb5,6\"]      <- -1\n",
    "a[\"Cdc20&Cdc14\",  \"Clb1,2\"]      <- -1\n",
    "a[\"Sic1\",         \"Clb1,2\"]      <- -1\n",
    "a[\"Cdh1\",         \"Clb1,2\"]      <- -1\n",
    "\n",
    "# define the autorepression edges (yellow loops)\n",
    "a[\"Cln3\",         \"Cln3\"]        <- -1\n",
    "a[\"Swi5\",         \"Swi5\"]        <- -1\n",
    "a[\"Cdc20&Cdc14\",  \"Cdc20&Cdc14\"] <- -1\n",
    "a[\"Cln1,2\",       \"Cln1,2\"]      <- -1\n",
    "a[\"Mcm1/SFF\",     \"Mcm1/SFF\"]    <- -1\n",
    "\n",
    "# define the fixed points of the model\n",
    "num.fp <- 7\n",
    "fixed.points <- matrix(nrow=num.fp,\n",
    "                       ncol=N)\n",
    "fixed.points[,] <- 0\n",
    "colnames(fixed.points) <- nodes\n",
    "\n",
    "# define the nonzero entries of the fixed-point matrix\n",
    "fixed.points[1,\"Cdh1\"] <- 1\n",
    "fixed.points[1,\"Sic1\"] <- 1\n",
    "fixed.points[2,\"SBF\"] <- 1\n",
    "fixed.points[2,\"Cln1,2\"] <- 1\n",
    "fixed.points[3,\"MBF\"] <- 1\n",
    "fixed.points[3,\"Cdh1\"] <- 1\n",
    "fixed.points[3,\"Sic1\"] <- 1\n",
    "fixed.points[4,\"Sic1\"] <- 1\n",
    "fixed.points[5,\"MBF\"] <- 1\n",
    "fixed.points[5,\"Sic1\"] <- 1\n",
    "fixed.points[7,\"Cdh1\"] <- 1\n",
    "\n",
    "# set up the vector for recording the number of times you end up at each\n",
    "# of the different fixed points\n",
    "basin.counts <- vector(length=num.fp, mode=\"integer\")\n",
    "basin.counts[] <- 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function `hamming.dist` that gives the hamming distance between two states of the Boolean network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "hamming.dist <- function(x1,x2) {\n",
    "# return the Hamming distance between state vectors x1 and x2\n",
    "  return(sum(abs(x1-x2)))  \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function `evolve` that takes the network from one Boolean vector state to another Boolean vector state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "evolve <- function(state) {\n",
    "  # \"state\" is a vector of length N, containing zeros or ones.\n",
    "  # A zero for a given node means it is in the \"off\" state; a\n",
    "  # one means that the node is in the \"on\" state. Update for one\n",
    "  # time-step according to the instructions in Eq. 1 of the Li et al. paper,\n",
    "                                        # and return the new state\n",
    "\n",
    "  # compute the influence function as product of state and the a matrix\n",
    "  result <- state %*% a\n",
    "\n",
    "  # positive entries get set to 1\n",
    "  result[result > 0] <- 1\n",
    "\n",
    "  # for any node for which the influence function is zero, use the current state\n",
    "  result[result == 0] <- state[result == 0]\n",
    "  \n",
    "  # negative entries get set to 0\n",
    "  result[result < 0] <- 0\n",
    "  \n",
    "  return(result)    \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write a function that runs 10,000 simulations of the network. In each simulation, the procedure is:\n",
    "- create a random binary vector of length 12, and call that vector `state` (make sure the 1st position of the vector is set to zero)\n",
    "- iteratively call \"evolve\", passing the `state` to evolve and then updating `state` with the return value from `evolve`\n",
    "- check if `state` changes in the last call to `evolve`; if it does not, then you have reached a fixed point; stop iterating\n",
    "- compare the state to the rows of `fixed.points`; for the unique row `i` for which you find a match, increment the element in position `i` of `basin.counts`\n",
    "- print out `basin.counts`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "basin_ids\n",
       "   1    2    3    4    5    6    7    8    9   10   11 \n",
       "4293  283 4536  395  197   21  180    1   59   17   18 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}