class20_booleannetwork_python3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CSX46 - Class 20 - Boolean Networks\n",
    "\n",
    "In this class, we will implement the *Li et al.* (PNAS, 2004) Boolean network model of the protein network controlling the yeast cell cycle. This model has 12 nodes and an interesting way of implementing state changes (dynamics) that differs slightly from a \"canonical\" Boolean network. After implementing the model, we will run a bunch of forward simulations using random initial states of the proteins, and we will determine (for each random initial state) which fixed-point (i.e., steady) state the system ends up in. We will compute the size of the \"basin\" for each fixed point (sometimes called the \"attractor basin\") as the number of random initial states that ultimately lead to that particular fixed-point state. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't need many packages for this notebook, just a few of the basics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import numpy\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will define a list of \"node names\" and an adjacency matrix. You will see that the ajacency matrix has some negative entries!  These correspond to inhibitory edges. I have saved you a lot of trouble by typing in the adjacency matrix; I made many mistakes before I got it right!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "nodes = ['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",
    "N = len(nodes)\n",
    "\n",
    "# define the transition matrix \n",
    "a = numpy.zeros([N, N])\n",
    "a[0,1] = 1\n",
    "a[1,1] = -1\n",
    "a[1,2] = 1\n",
    "a[1,11] = 1\n",
    "a[2,3] = 1\n",
    "a[3,4] = 1\n",
    "a[3,6] = -1\n",
    "a[3,7] = 1\n",
    "a[3,9] = -1\n",
    "a[4,4] = -1\n",
    "a[4,5] = 1\n",
    "a[4,7] = 1\n",
    "a[4,8] = 1\n",
    "a[5,5] = -1\n",
    "a[5,6] = 1\n",
    "a[6,3] = -1\n",
    "a[6,7] = -1\n",
    "a[7,2] = -1\n",
    "a[7,4] = 1\n",
    "a[7,5] = -1\n",
    "a[7,6] = -1\n",
    "a[7,8] = 1\n",
    "a[7,9] = -1\n",
    "a[7,11] = -1\n",
    "a[8,3] = -1\n",
    "a[8,5] = 1\n",
    "a[8,6] = 1\n",
    "a[8,7] = -1\n",
    "a[8,8] = -1\n",
    "a[8,9] = 1\n",
    "a[9,7] = -1\n",
    "a[10,6] = -1\n",
    "a[10,9] = -1\n",
    "a[10,10] = -1\n",
    "a[11,10] = 1\n",
    "a = numpy.matrix(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "        [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],\n",
       "        [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "        [ 0.,  0.,  0.,  0.,  1.,  0., -1.,  1.,  0., -1.,  0.,  0.],\n",
       "        [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  1.,  1.,  0.,  0.,  0.],\n",
       "        [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],\n",
       "        [ 0.,  0.,  0., -1.,  0.,  0.,  0., -1.,  0.,  0.,  0.,  0.],\n",
       "        [ 0.,  0., -1.,  0.,  1., -1., -1.,  0.,  1., -1.,  0., -1.],\n",
       "        [ 0.,  0.,  0., -1.,  0.,  1.,  1., -1., -1.,  1.,  0.,  0.],\n",
       "        [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  0.,  0.,  0.,  0.],\n",
       "        [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  0.,  0., -1., -1.,  0.],\n",
       "        [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Li *et al.* article describes seven \"fixed points\" for their Boolean network model. They can be represented as a 7x12 matrix, as shown here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the matrix of states for the fixed points\n",
    "num_fp = 7\n",
    "fixed_points = numpy.zeros([num_fp, N])\n",
    "fixed_points[0, 6] = 1\n",
    "fixed_points[0, 9] = 1\n",
    "fixed_points[1, 10] = 1\n",
    "fixed_points[1, 11] = 1\n",
    "fixed_points[2, 2] = 1\n",
    "fixed_points[2, 6] = 1\n",
    "fixed_points[2, 9] = 1\n",
    "fixed_points[3, 6] = 1\n",
    "fixed_points[4, 2] = 1\n",
    "fixed_points[4, 6] = 1\n",
    "fixed_points[6, 9] = 1\n",
    "fixed_points = numpy.matrix(fixed_points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the fixed-point matrix; each row is a Boolean \"state vector\" for the 12 nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0., 0., 0., 0., 0., 0., 1., 0., 0., 1., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1.],\n",
       "        [0., 0., 1., 0., 0., 0., 1., 0., 0., 1., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fixed_points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function `hamming.dist` that gives the hamming distance between two states of the Boolean network (as numpy arrays of ones and zeroes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hamming_dist(x1, x2):\n",
    "    return numpy.sum(numpy.abs(x1-x2))"
   ]
  },
  {
   "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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evolve(state):\n",
    "    result = numpy.array(a.transpose().dot(state))\n",
    "    result = numpy.reshape(result, N)\n",
    "    result[result > 0] = 1\n",
    "    result[result == 0] = state[result == 0]\n",
    "    result[result < 0] = 0\n",
    "    return result"
   ]
  },
  {
   "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 zeroth element 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 `j` for which you find a match, add `j` to the list `basin_ids`.\n",
    "- use `numpy.bincount` to count the number of times each basin ID appears in `basin_ids`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "basin_ids = []\n",
    "for _ in itertools.repeat(None, 10000):\n",
    "    state = [0]\n",
    "    for pos in range(0, (N-1)):\n",
    "        state.append(random.randint(0,1))\n",
    "    state = numpy.array(state)\n",
    "    state_new = numpy.array([-1]*N)\n",
    "    while(True):\n",
    "        state_new = evolve(state)\n",
    "        if hamming_dist(state, state_new) == 0:\n",
    "            break\n",
    "        state = state_new\n",
    "    for j in range(0, num_fp):\n",
    "        fp_state = numpy.array(fixed_points[j,])\n",
    "        fp_state = numpy.reshape(fp_state, N)\n",
    "        if hamming_dist(state, fp_state) == 0:\n",
    "            basin_ids.append(j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([8638,  714,  527,   41,   34,   37,    9])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fp_counts = numpy.bincount(basin_ids)\n",
    "fp_counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which fixed-point has the highest count?  Compare to Table 1 in Li *et al.* 2004."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Cell Size': 0.0,\n",
       " 'Cln3': 0.0,\n",
       " 'MBF': 0.0,\n",
       " 'Clb5,6': 0.0,\n",
       " 'Mcm1/SFF': 0.0,\n",
       " 'Swi5': 0.0,\n",
       " 'Sic1': 1.0,\n",
       " 'Clb1,2': 0.0,\n",
       " 'Cdc20&Cdc14': 0.0,\n",
       " 'Cdh1': 1.0,\n",
       " 'Cln1,2': 0.0,\n",
       " 'SBF': 0.0}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "{nodes[i]: fixed_points[numpy.argmax(fp_counts),].tolist()[0][i] for i in range(0,12)}"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}