{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from random import random\n",
    "\n",
    "L = 60\n",
    "\n",
    "def rand_i ():\n",
    "    return int(random() * L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start with simulating the basic 3-state nucleosome model for a system with N=60 nucleosomes.\n",
    "\n",
    "Then do variants Q7.2.1, Q7.2.2.\n",
    "\n",
    "First where U-->M and U-->A local (only acting on nearest neighbors), whereas recruitments A--> U and M-->U are global (act randomly in entire system). This local-global model variant needs two parameters:\n",
    "the parameter alpha for how often one chose recruitment (should be around 90 for random conversions), and a parameter beta=0.5 for the probability that a given recruitment attempt should be global. When attempting a global attempt, select one recruiting nucleosome and another that it acts on anywhere in the system. The move can only be performed if the two chosen nuleosomes are of type M and A, or A and M. Similarly a local move involves selection of two neighbor nucleosomes where the nucleosome to be changed have to start in U state. Use again an L=60 system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 42956 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "steps = 100000\n",
    "alpha = .9\n",
    "\n",
    "nuc = np.zeros(L)\n",
    "nucs_A = np.zeros(steps)\n",
    "nucs_M = np.zeros(steps)\n",
    "\n",
    "for n in range(steps):\n",
    "    i = rand_i()\n",
    "\n",
    "    if (random() <= alpha):\n",
    "        # Recruited conversion\n",
    "        j = i\n",
    "        while j == i:\n",
    "            j = rand_i()\n",
    "        \n",
    "        nuc[i] += nuc[j]\n",
    "    \n",
    "    else:\n",
    "        # Noisy change\n",
    "        if nuc[i] == 0:\n",
    "            if random() < 0.5:\n",
    "                nuc[i] = 1\n",
    "            else:\n",
    "                nuc[i] = -1\n",
    "        else:\n",
    "            nuc[i] = 0\n",
    "    \n",
    "    if nuc[i]:\n",
    "        nuc[i] = nuc[i] / np.abs(nuc[i])\n",
    "    \n",
    "    nucs_A[n] = np.count_nonzero(nuc == -1)\n",
    "    nucs_M[n] = np.count_nonzero(nuc == 1)\n",
    "\n",
    "\n",
    "plt.title(\"Global\")\n",
    "plt.plot(nucs_A, label=\"A\")\n",
    "plt.plot(nucs_M, label=\"M\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 68712 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "steps = 100000\n",
    "alpha = .9\n",
    "beta = .5\n",
    "\n",
    "nuc = np.zeros(L)\n",
    "nucs_A = np.zeros(steps)\n",
    "nucs_M = np.zeros(steps)\n",
    "\n",
    "for n in range(steps):\n",
    "    i = rand_i()\n",
    "\n",
    "    if random() <= alpha:\n",
    "        # Recruited conversion\n",
    "\n",
    "        if random() <= beta:\n",
    "            # global\n",
    "            j = i\n",
    "            while j == i:\n",
    "                j = rand_i()\n",
    "            \n",
    "            # only A->U, M->U\n",
    "            if nuc[i] + nuc[j] == 0:\n",
    "                nuc[i] = 0\n",
    "        else:\n",
    "            # local\n",
    "            if i == 0:\n",
    "                j = 1\n",
    "            elif i == L-1:\n",
    "                j = L-2\n",
    "            elif random() < .5:\n",
    "                j = i+1\n",
    "            else:\n",
    "                j = i-1\n",
    "\n",
    "            # only U->A, U->M\n",
    "            if np.abs(nuc[i] + nuc[j]) == 1:\n",
    "                nuc[i] += nuc[j]\n",
    "\n",
    "    \n",
    "    else:\n",
    "        # Noisy change\n",
    "        if nuc[i] == 0:\n",
    "            if random() < 0.5:\n",
    "                nuc[i] = 1\n",
    "            else:\n",
    "                nuc[i] = -1\n",
    "        else:\n",
    "            nuc[i] = 0\n",
    "    \n",
    "    if nuc[i]:\n",
    "        nuc[i] = nuc[i] / np.abs(nuc[i])\n",
    "    \n",
    "    nucs_A[n] = np.count_nonzero(nuc == -1)\n",
    "    nucs_M[n] = np.count_nonzero(nuc == 1)\n",
    "\n",
    "\n",
    "plt.title(\"A->U, M->U global and U->A, U->M local\")\n",
    "plt.plot(nucs_A, label=\"A\")\n",
    "plt.plot(nucs_M, label=\"M\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 63816 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "steps = 100000\n",
    "alpha = .9\n",
    "beta = .5\n",
    "\n",
    "nuc = np.zeros(L)\n",
    "nucs_A = np.zeros(steps)\n",
    "nucs_M = np.zeros(steps)\n",
    "\n",
    "for n in range(steps):\n",
    "    i = rand_i()\n",
    "\n",
    "    if random() <= alpha:\n",
    "        # Recruited conversion\n",
    "\n",
    "        if random() <= beta:\n",
    "            # global\n",
    "            j = i\n",
    "            while j == i:\n",
    "                j = rand_i()\n",
    "            \n",
    "            # only U->A, U->M\n",
    "            if np.abs(nuc[i] + nuc[j]) == 1:\n",
    "                nuc[i] += nuc[j]\n",
    "           \n",
    "            \n",
    "        else:\n",
    "            # local\n",
    "            if i == 0:\n",
    "                j = 1\n",
    "            elif i == L-1:\n",
    "                j = L-2\n",
    "            elif random() < .5:\n",
    "                j = i+1\n",
    "            else:\n",
    "                j = i-1\n",
    "            \n",
    "            # only A->U, M->U\n",
    "            if nuc[i] + nuc[j] == 0:\n",
    "                nuc[i] = 0\n",
    "\n",
    "    \n",
    "    else:\n",
    "        # Noisy change\n",
    "        if nuc[i] == 0:\n",
    "            if random() < 0.5:\n",
    "                nuc[i] = 1\n",
    "            else:\n",
    "                nuc[i] = -1\n",
    "        else:\n",
    "            nuc[i] = 0\n",
    "    \n",
    "    if nuc[i]:\n",
    "        nuc[i] = nuc[i] / np.abs(nuc[i])\n",
    "    \n",
    "    nucs_A[n] = np.count_nonzero(nuc == -1)\n",
    "    nucs_M[n] = np.count_nonzero(nuc == 1)\n",
    "\n",
    "\n",
    "plt.title(\"U->A, U->M global and A->U, M->U local\")\n",
    "plt.plot(nucs_A, label=\"A\")\n",
    "plt.plot(nucs_M, label=\"M\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
