{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from math import sqrt\n",
    "\n",
    "def pq (a,b,c):\n",
    "    return (-b + sqrt(b**2 - 4*a*c)) / 2*a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-3.792751279232142\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 47356 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Transcription regulation\n",
    "# Question 3.2.1, 3.2.2\n",
    "\n",
    "dt = 0.01\n",
    "t = np.linspace(0,20,int(20 / dt))\n",
    "C1 = np.zeros_like(t)\n",
    "C2 = np.zeros_like(t)\n",
    "C3 = np.zeros_like(t)\n",
    "C4 = np.zeros_like(t)\n",
    "m = np.zeros_like(t)\n",
    "C5 = np.zeros_like(t)\n",
    "tau = 10\n",
    "\n",
    "R = np.where(t <= 10, .1, 10)\n",
    "\n",
    "C1[0] = tau / (R[0] + 1)\n",
    "C2[0] = 1 / (R[0] + 1 / tau)\n",
    "C3[0] = tau * R[0] / (R[0] + 1)\n",
    "C4[0] = pq(10,1,-2 * tau * R[0]/(R[0] + 1))\n",
    "\n",
    "print((R[0] / (R[0] + 1) * 2 / (10 * C4[0] + 1) - C4[0] / tau))\n",
    "\n",
    "m[0] = 1 / (R[0] + 1 / tau)\n",
    "C5[0] = R[0] * m[0]\n",
    "\n",
    "for i in range (0, int(20 / dt) -1):\n",
    "    C1[i+1] = C1[i] + (1 / (1 + R[i]) - C1[i] / tau)  * dt\n",
    "    C2[i+1] = C2[i] + (1 - (R[i] + 1  / tau) * C2[i])  * dt\n",
    "    C3[i+1] = C3[i] + (R[i] / (R[i] + 1) - C3[i] / tau)  * dt\n",
    "    C4[i+1] = C4[i] + (R[i] / (R[i] + 1) * 2 / (10 * C4[i] + 1) - C4[i] / tau) * dt\n",
    "    m[i+1] = m[i] + (1 - m[i] / tau - R[i] * m[i]) * dt\n",
    "    C5[i+1] = C5[i] + (R[i] * m[i] - C5[i] ) * dt\n",
    "\n",
    "\n",
    "plt.title(\"R increases by factor 100\")\n",
    "plt.plot(t, C1, label=\"transcriptional repression\")\n",
    "plt.plot(t, C2, label=\"protein degradation\")\n",
    "plt.plot(t, C3, label=\"transcriptional activation\")\n",
    "#plt.plot(t, C4, label=\"self repression\")\n",
    "plt.plot(t, C5, label=\"activated reservoir\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fd56e430ee0>]"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 19932 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# sRNA regulation\n",
    "# Questions 3.3.1, 3.3.2, 3.3.3, 3.3.4\n",
    "\n",
    "dt = 0.01\n",
    "t = np.linspace(0,20,int(20 / dt))\n",
    "\n",
    "m = np.zeros_like(t)\n",
    "s = np.zeros_like(t)\n",
    "tau = 1\n",
    "\n",
    "gamma = 10\n",
    "alpha = np.where(t <= 10, .1, 4)\n",
    "\n",
    "\n",
    "m[0] = pq(- gamma / tau, (gamma - alpha[0] * gamma - 1 / tau), 1)\n",
    "s[0] = alpha[0] / (1 + gamma * m[0])\n",
    "\n",
    "\n",
    "for i in range (0, int(20 / dt) -1):\n",
    "    m[i+1] = m[i] + (1 - gamma * alpha[i] * s[i] - m[i] / tau) * dt\n",
    "    s[i+1] = s[i] + (alpha[i] - gamma * alpha[i] * s[i] - s[i]) * dt\n",
    "\n",
    "plt.plot(t, m)\n",
    "plt.plot(t, s)\n",
    "\n",
    "# FIXME\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
}
