biophys_calc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Self-Assembly of a Dimer System -- Companion Notebook**\n",
    "\n",
    "# Biophysical Parameter Calculations "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import *\n",
    "import pandas as pd\n",
    "from scipy.special import lambertw\n",
    "from decimal import Decimal\n",
    "from random import randint\n",
    "import time\n",
    "\n",
    "start_clock = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Contents\n",
    "\n",
    "**0.** Calculating $E_0$ and $\\Delta$ from binding free energies\n",
    "\n",
    "**I.** ssDNA - ssDNA interaction calculations\n",
    "\n",
    "- (a) Free Energy Estimate\n",
    "\n",
    "- (b) Calculating $(N)_{\\text{max}}$, $(NV)_{\\text{max}}$, and Type I inequality\n",
    "\n",
    "**II.** TF-DNA Free interaction calculations\n",
    "\n",
    "- (a) Free Energy Estimate\n",
    "\n",
    "- (b) Calculating $(N)_{\\text{max}}$, $(NV)_{\\text{max}}$, and Type I inequality\n",
    "\n",
    "**III.** Protein-protein interaction calculations\n",
    "\n",
    "- (a) Free Energy Estimate\n",
    "\n",
    "- (b) Calculating $(N)_{\\text{max}}$, $(NV)_{\\text{max}}$, and Type I inequality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0. Calculating $E_0$ and $\\Delta$ from binding free energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the subsequent code, we determine the energy parameters $E_0$ and $\\Delta$ from estimates of biomolecular binding free energies as follows. \n",
    "\n",
    "\n",
    "For a binding/association reaction, we define changes in Gibbs free energy as $\\delta G$ and changes in entropy as $\\delta S$. Thermodynamically, Gibbs free energy is distinct from Helmholtz free energy the latter of which is defined in terms of the former as $F = G + PV$, where $p$ and $V$ are system pressure and volume respectively. However, for systems at constant pressure and volume (as is presumed for our biomolecular dimerization systems), $\\delta p$ and $\\delta V$ are both zero, and changes in Helmholtz free energy equal changes in Gibbs free energy, $\\delta G = \\delta F$. Therefore we take $\\delta G = \\delta E - T \\delta S$, where $\\delta E$ is the change in energy upon binding. Solving for $\\delta E$, we then obtain \n",
    "\n",
    "$$\\delta E = \\delta G + T\\delta S.$$\n",
    "\n",
    "In particular, for the model, the change in energy for a correct/specific interaction is $(\\delta E)_{\\text{correct}} = -(E_0+\\Delta)$ and the change in energy for an incorrect/nonspecific interaction is $(\\delta E)_{\\text{incorrect}} = - E_0$. Therefore, we have \n",
    "\n",
    "\\begin{eqnarray}\n",
    "-(E_0+\\Delta) & = (\\delta G)_{\\text{correct}} + T(\\delta S)_{\\text{incorrect}},\\\\\n",
    "-E_0 & = (\\delta G)_{\\text{incorrect}} + T(\\delta S)_{\\text{incorrect}}.\n",
    "\\end{eqnarray}\n",
    "\n",
    "In our model the only contribution to the change in entropy upon binding of two monomers is the change in translational entropy, and this change in translational entropy is independent of whether the binding is correct or incorrect. Therefore, we take \n",
    "\n",
    "$$ (\\delta S)_{\\text{correct}}  = (\\delta S)_{\\text{incorrect}} \\equiv (\\delta S)_{\\text{trans}} $$\n",
    "\n",
    "and the equations become \n",
    "\n",
    "\\begin{eqnarray}\n",
    "-(E_0+\\Delta) & = (\\delta G)_{\\text{correct}} + T(\\delta S)_{\\text{trans}}, \\\\\n",
    "-E_0 & = (\\delta G)_{\\text{incorrect}} + T(\\delta S)_{\\text{trans}}.\n",
    "\\end{eqnarray}\n",
    "\n",
    "From these two equations, we find that $\\Delta$ can be computed from the equation \n",
    "\n",
    "$$ \\Delta = (\\delta G)_{\\text{incorrect}}- (\\delta G)_{\\text{correct}},$$\n",
    "\n",
    "and $E_0$ is computed from \n",
    "\n",
    "$$ E_0 =  -(\\delta G)_{\\text{incorrect}} - T(\\delta S)_{\\text{trans}}.$$\n",
    "\n",
    "The binding free energy $(\\delta G)_{\\text{incorrect}}$ is what we obtain from the biophysical data. We obtain $T(\\delta S)_{\\text{trans}}$ from the corresponding ideal gas expression. Specifically, the Sackur-Tetrode equation states that the entropy of an ideal solution (written in kcal/mol) is\n",
    "\n",
    "\\begin{equation}\n",
    "T S/N_A = R T \\ln \\left[ e^{5/2}\\frac{V_{\\text{st}}}{N_A\\lambda^3} \\right]\n",
    "\\end{equation}\n",
    "\n",
    "where $\\lambda = h/\\sqrt{2\\pi m_0 k_BT}$ is the thermal de Broglie wavelength of the particle of interest, $N_A = 6.022 \\times 10^{23}$ Avogadro's number,  $V_{\\text{st}} = $ 1 L (standard volume of ideal solutions), and $R = 1.987\\times10^{-3}$ kcal/K$^{-1}$mol$^{-1}$. We define the entropy in terms of $V_{\\text{st}}$ and $N_A$ in order to match the conditions under which the binding free energy is defined; $\\delta G$ is defined in units of kcal/mol i.e., in units of kcal per Avogadro's number of particles in a 1 L volume. \n",
    "\n",
    "As a general example, consider two motile monomers $A$ and $B$ which form the dimer $AB$. If the corresponding thermal de Broglie wavelengths of $A$, $B$, and $AB$ are $\\lambda_{\\alpha}$, $\\lambda_{\\beta}$ and $\\lambda_{\\alpha\\beta}$, then the change in entropy upon binding is  \n",
    "\n",
    "\\begin{eqnarray}\n",
    "T \\delta S_{\\text{trans}}/ N_{A} = T (S_{\\alpha\\beta}- S_{\\alpha} - S_{\\beta})/N_A = -RT \\ln \\left[e^{5/2}\\frac{V_{\\text{st}}}{\\lambda_{\\alpha}^3}\\frac{\\lambda_{\\alpha\\beta}^3}{\\lambda_{\\beta}^3}  \\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "We later apply this formula to calculate the binding energy parameters for our three systems of interest. \n",
    "\n",
    "It is known that for real protein systems, the ideal-solution approximation for entropy leads to an overestimate in the change in translational entropy upon binding since real protein monomers rarely explore their entire surrounding volume [1-2]. Thus, we can interpret our chosen entropy as providing an upper limit on the translational entropy contribution to the free energy and thus a lower limit on the obtained binding energy parameters in the model. \n",
    "\n",
    "\n",
    "\n",
    "- [1] Murphy, K. P., Xie, D., Thompson, K. S., Amzel, L. M., & Freire, E. (1994). Entropy in biological binding processes: estimation of translational entropy loss. *Proteins: Structure, Function, and Bioinformatics*, 18(1), 63-67.\n",
    "\n",
    "- [2] Amzel, L. M. (1997). Loss of translational entropy in binding, folding, and catalysis. *Proteins: Structure, Function, and Bioinformatics*, 28(2), 144-149.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___________\n",
    "___________"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### I. ssDNA-ssDNA interaction calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (a). Free Energy Estimate\n",
    "\n",
    "We develop code to calculate the free energy of binding for ssDNA+ssDNA(complement) to dsDNA given the sequence of one of the strands. \n",
    "\n",
    "This code is based on the results in \n",
    "\n",
    "* SantaLucia, J. (1998). A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics. *Proceedings of the National Academy of Sciences, 95(4), 1460-1465.\n",
    "\n",
    "The main equation we will be using is \n",
    "\n",
    "$$ \\Delta G^{\\circ}(\\text{total}) = \\sum_{i}n_i \\Delta G^{\\circ}(i) + \\Delta G^{\\circ}(\\text{with term./ A-T})+\\Delta G^{\\circ}(\\text{with term./ G-C})+ \\Delta G^{\\circ}(\\text{sym.})$$\n",
    "\n",
    "where $n_i$ is the number of nearest-neighbor base-pair interactions of type $i$, and $\\Delta G^{\\circ}(i)$ is the corresponding free energy; $\\Delta G^{\\circ}(\\text{with term./ A-T})$ and $\\Delta G^{\\circ}(\\text{with term./ G-C})$ are the additional free energies added for double strands ending in A-T and G-C respectively; $\\Delta G^{\\circ}(\\text{sym.})$ is added to the result for self-complementary sequences. This equation applies for a 1M Na$^+$ solution.\n",
    "\n",
    "Here is the table of free energies for these interactions at $37^{\\circ}$ C. \n",
    "\n",
    "| Sequence | $\\Delta G(i)$, kcal/mol|\n",
    "| :---: | :---:  |\n",
    "| AA/TT | -1.00  |\n",
    "| AT/TA | -0.88  |\n",
    "| TA/AT | -0.58  |\n",
    "| CA/GT | -1.45  |\n",
    "| GT/CA | -1.44  |\n",
    "| CT/GA | -1.28  |\n",
    "| GA/CT | -1.30  |\n",
    "| CG/GC | -2.17  |\n",
    "| GC/CG | -2.24  |\n",
    "| GG/CC | -1.84  |\n",
    "|Terminal G-C| 0.98  | \n",
    "|Terminal A-T| 1.03  |\n",
    "\n",
    "**NOTE:** For a system with an arbitrary salt concentration we use the formula \n",
    "\n",
    "$$ \\Delta G^{\\circ}_{37}(\\text{oligomer}, [\\text{Na}^{+}]) = \\Delta G^{\\circ}_{37}(\\text{oligomer}, 1 \\text{ M NaCl}) - 0.114 \\times N \\times \\ln [ \\text{Na}^{+}], $$ \n",
    "\n",
    "where $N$ is the number of total phosphates in the duplex divided by 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Preliminary definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Takes a 5' to 3' sequence and outputs \n",
    "# the complementary sequence in the 5' to 3' direction\n",
    "\n",
    "def reverse(s): \n",
    "    st = \"\" \n",
    "    for i in s: \n",
    "        st = i + st\n",
    "    return st\n",
    "\n",
    "def comp_seq(string):\n",
    "    \n",
    "    empt_str = ''\n",
    "    for k in string:\n",
    "        if k == 'A':\n",
    "            \n",
    "            empt_str = empt_str +'T'\n",
    "        elif k == 'T':\n",
    "            \n",
    "            empt_str = empt_str + 'A'\n",
    "        \n",
    "        elif k == 'G':\n",
    "            \n",
    "            empt_str = empt_str +'C'\n",
    "            \n",
    "        elif k == 'C':\n",
    "            \n",
    "            empt_str = empt_str +'G'\n",
    "            \n",
    "     \n",
    "    return(reverse(empt_str))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## dictionary of energies for various sequences\n",
    "## From Last column of Table 1 of referenced paper\n",
    "\n",
    "en_dict = {\n",
    "    \"AA\": -1.00,\n",
    "    \"TT\": -1.00,\n",
    "    \"AT\": -0.88,\n",
    "    \"TA\": -0.58,\n",
    "    \"AG\": -1.28,\n",
    "    \"AC\": -1.44,\n",
    "    \"TG\": -1.45,\n",
    "    \"TC\": -1.30,\n",
    "    \"CA\": -1.45,\n",
    "    \"GT\": -1.44,\n",
    "    \"CT\": -1.28,\n",
    "    \"GA\": -1.30,\n",
    "    \"CG\": -2.17,\n",
    "    \"GC\": -2.24,\n",
    "    \"GG\": -1.84,\n",
    "    \"CC\": -1.84\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.24"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# testing dictionary\n",
    "en_dict['GC']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### ΔG of binding in 1M salt solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "## computes free energy in kcal/mol of DNA duplex in a 1M salt solution\n",
    "## From Eq. 1 of referenced paper\n",
    "\n",
    "\n",
    "def delG(string):\n",
    "    \n",
    "    Gtot = 0\n",
    "    \n",
    "    init = string[0]\n",
    "    fin = string[len(string)-1]\n",
    "    \n",
    "    for k in range(len(string)-1):\n",
    "\n",
    "        min_str = string[k:k+2]\n",
    "        Gtot = Gtot + en_dict[min_str]\n",
    "    \n",
    "    # A-T terminal pair additional energy\n",
    "    if fin == 'A' or fin == 'T':\n",
    "        \n",
    "        Gtot = Gtot + 1.03\n",
    "        \n",
    "    # G-C terminal pair additional energy\n",
    "    elif fin == 'G' or fin == 'C':\n",
    "        \n",
    "        Gtot = Gtot + 0.98\n",
    "        \n",
    "    # A-T initial pair additional energy\n",
    "    if init == 'A' or init == 'T':\n",
    "        \n",
    "        Gtot = Gtot + 1.03\n",
    "        \n",
    "    \n",
    "    # G-C initial pair additional energy\n",
    "    elif init == 'G' or init == 'C':\n",
    "        \n",
    "        Gtot = Gtot + 0.98\n",
    "        \n",
    "    ## additional free energy for self-complementarity\n",
    "    if string == reverse(string):\n",
    "        \n",
    "        Gtot = Gtot + 0.43\n",
    "        \n",
    "    return Gtot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-5.35"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## testing calculation on sequence and corresponding value found in paper \n",
    "delG('CGTTGA')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-7.105427357601002e-15"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## testing that delta G gives same energy for sequence and complement\n",
    "\n",
    "test_seq = 'CGCTAGCATCAGTCAGTCAGT'\n",
    "delG(test_seq)- delG(comp_seq(test_seq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "## randomly generate 100000 20-base-pair sequences for ssDNA\n",
    "\n",
    "bp_list = 'CGTA'\n",
    "\n",
    "leng = 100000\n",
    "seq_vec = ['']*leng\n",
    "for k in range(leng):\n",
    "    \n",
    "    for j in range(20):\n",
    "        u = randint(0,3)\n",
    "        \n",
    "        seq_vec[k] = seq_vec[k]+bp_list[u]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-26.520000000000007"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test; compare with result in http://biophysics.idtdna.com/\n",
    "delG('CGCTAGCATCAGTCAGTCAGT')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "## computing delta G for our list of sequences \n",
    "## in 1M salt solution\n",
    "\n",
    "delG_list = np.ones(leng)\n",
    "\n",
    "for k in range(leng):\n",
    "    \n",
    "    delG_list[k] = delG(seq_vec[k])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean ΔGs (1M salt): -24.692644100000003\n",
      "Standard Deviation of ΔGs (1M salt): 2.5669832418103544\n"
     ]
    }
   ],
   "source": [
    "## computing mean and variance of ΔGs in 1M salt solution\n",
    "\n",
    "print('Mean ΔGs (1M salt):',np.mean(delG_list))\n",
    "print('Standard Deviation of ΔGs (1M salt):', np.std(delG_list))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### ΔG of binding in arbitary salt solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "## calculation of ΔG with an incoporation of salt concentration \n",
    "## From Eq. 7 of referenced paper\n",
    "\n",
    "def delG_salt(string, salt_concent):\n",
    "    \n",
    "    # salt_concentration is in units of molarity\n",
    "    # e.g., 1M -> salt_concent  = 1\n",
    "    \n",
    "    # number of phosphates in the duplex divided by 2\n",
    "    # assumes each nucleotide has a phosphate\n",
    "    N = len(string)\n",
    "    \n",
    "    return delG(string) -0.114*N*np.log(salt_concent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-19.348216937111744"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test; compare with result in http://biophysics.idtdna.com/\n",
    "delG_salt('ACTGACTGACTGATGCTAGCG', 0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "## computing delta G for our list of sequences \n",
    "## in 1M salt solution\n",
    "\n",
    "delG_list_salt = np.ones(leng)\n",
    "\n",
    "for k in range(leng):\n",
    "    \n",
    "    delG_list_salt[k] = delG_salt(seq_vec[k], 0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean ΔGs (50 mM salt): -17.8623745162969\n",
      "Standard Deviation of ΔGs (50 mM salt): 2.5669832418103544\n"
     ]
    }
   ],
   "source": [
    "## computing mean and variance of ΔGs in 50 mM salt solution\n",
    "\n",
    "print('Mean ΔGs (50 mM salt):',np.mean(delG_list_salt))\n",
    "print('Standard Deviation of ΔGs (50 mM salt):', np.std(delG_list_salt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (b). Calculating $(N)_{\\text{max}}$, $(NV)_{\\text{max}}$, and Type I inequality\n",
    "\n",
    "We calculate $(N)_{\\text{max}}$, $(NV)_{\\text{max}}$, and the right-hand-side of the \"Type I inequality\" for the TF-DNA system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Relevant Inequalities**\n",
    "\n",
    "$$\n",
    "(NV)_{\\text{max}} = \\sqrt{2} \\lambda_{0}^3 e^{\\beta(\\Delta+E_0)}, \\qquad \\text{[$T= T_{\\text{I}}$]}\n",
    "$$\n",
    "where $\\lambda_{0} = h/\\sqrt{2\\pi m_0 k_BT}$\n",
    "\n",
    "\\begin{equation}\n",
    "(N)_{\\text{max}} = e^{\\beta \\Delta}/2, \\qquad \\text{[$T= T_{\\text{II}}$]} \n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "2N < \\exp\\left[\\frac{3\\Delta}{2E_0} W_{0}\\left(\\frac{E_0}{3E_V}\\right)\\right], \\qquad \\text{[Type I inequality, $T_{\\text{I}}< T_{\\text{II}}$]}\n",
    "\\end{equation}\n",
    "\n",
    "where $E_{ V} = h^2/2\\pi m_0 V^{2/3}$.\n",
    "\n",
    "**Change in Translational Entropy upon Binding**\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "T\\left(S_{\\text{dsDNA}} - 2\\,S_{\\text{ssDNA}}\\right)/N_A & = R T \\ln \\left[e^{5/2}\\frac{V_{\\text{st}}}{N_A(\\lambda_0/\\sqrt{2})^3} \\right] - 2 R T \\ln \\left[e^{5/2}\\frac{V_{\\text{st}}}{N_A\\lambda_0^3} \\right]\\\\\n",
    "& = -R T \\ln \\left[e^{5/2}\\frac{V_{\\text{st}}}{2\\sqrt{2}\\,N_A\\lambda_0^3} \\right] \\\\\n",
    "& = -13.6 \\text{ kcal/mol}.\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-13.601045435350864"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Computing translational entropy contribution to free energy\n",
    "\n",
    "# We use SI units throughout, except for energies which are defined in kcal/mol\n",
    "\n",
    "Tphys = 310.15 # physiological temperature\n",
    "RTphys = 0.001987*Tphys # Gas constant times temperature\n",
    "V_stan = 1e-3 # 1 liter in meters cubed\n",
    "kB = 1.38*10**(-23) # Boltzmann's constant\n",
    "mass = 1.0793504*10**(-23) # 6500 Daltons mass of 20-base strand of ssDNA\n",
    "Na = 6.022*10**23 # Avogadro's number\n",
    "h = 6.626*10**(-34) # Planck's constant\n",
    "lam_0 = h/np.sqrt(2*np.pi*mass*kB*Tphys) # thermal deBroglie wavelength\n",
    "\n",
    "-RTphys*np.log(np.exp(5/2)*V_stan/(2*np.sqrt(2)*Na*lam_0**3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Energy Parameter Estimation**\n",
    "\n",
    "From the specific free energy of binding we find: \n",
    "\n",
    "$$ \\delta E = \\delta G + T \\delta S = -17.9 \\text{ kcal/mol} - 13.6 \\text{ kcal/mol} = -31.5  \\text{ kcal/mol} = -(E_0+\\Delta) =  -\\Delta$$\n",
    "\n",
    "where we took $E_0= 0$ in the final equality. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(N)max = 7.90E+21\n",
      "(NV)max= 4.16E+04\n"
     ]
    }
   ],
   "source": [
    "## Calculation (N)max, (NV)max, and Type I inequality\n",
    "\n",
    "# System paramaters for DNA system\n",
    "kB = 1.38*10**(-23) # boltzmann's constant\n",
    "mass = 1.0793504*10**(-23) # 6500 Daltons mass of 20-base strand of ssDNA\n",
    "Na = 6.022*10**23 # Avogadro's number\n",
    "h = 6.626*10**(-34) # Planck's constant\n",
    "\n",
    "E0 = 1e-5 # E0=0, \n",
    "Del = 31.5\n",
    "Tphys = 310.15 # physiological temperature\n",
    "RTphys = 0.001987*Tphys # Gas constant times temperature\n",
    "Vol = 1e-18 # we assume a volume of (1 micrometer)^3\n",
    "\n",
    "lam_0_phys = h/np.sqrt(2*np.pi*mass*kB*Tphys) # thermal deBroglie wavelength at phys temperature\n",
    "\n",
    "Ev = h**2/(2*np.pi*mass*Vol**(2/3))*(2.388e-4)*Na # Ev written in kcal/mol\n",
    "\n",
    "Nmax = np.exp(Del/RTphys)/2\n",
    "NVmax = np.sqrt(2)*lam_0_phys**3*np.exp((E0+Del)/RTphys)/(1e-18)\n",
    "# Type1_ineq = np.exp((3*Del)/(2*E0)*lambertw(E0/(3*Ev))).real  #overflow error\n",
    "\n",
    "print(\"(N)max =\", '%.2E' % Decimal(Nmax))\n",
    "print(\"(NV)max=\", '%.2E' % Decimal(NVmax))\n",
    "#print(\"exp((3Δ/2E0) W0(E0/3EV))=\", '%.2E' % Decimal(Type1_ineq))#overflow error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'3.38E+13'"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Argument of Type1_ineq exponential\n",
    "'%.2E' % Decimal(Del/(h**2/(2*np.pi*mass*Vol**(2/3))*(2.388e-4)*Na))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### I. TF-DNA interaction calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (a). Free Energy Estimate\n",
    "\n",
    "We estimate the specific and non-specific free energy of binding between TFs and DNA for some representative interactions\n",
    "\n",
    "Data is taken from Table 1 in \n",
    "\n",
    "* Jen‐Jacobson, L. (1997). Protein—DNA recognition complexes: Conservation of structure and binding energy in the transition state. *Biopolymers: Original Research on Biomolecules*, 44(2), 153-180."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 12 TFs and their specific and non-specific binding energies to DNA\n",
    "'''\n",
    "Transcribed from Table 1 in:\n",
    "\n",
    "Jen‐Jacobson, Linda. \n",
    "\"Protein—DNA recognition complexes: Conservation of structure \n",
    "and binding energy in the transition state.\" \n",
    "Biopolymers: Original Research on Biomolecules 44.2 (1997): 153-180.\n",
    "'''\n",
    "protein_DNA_dfs = pd.read_csv('data_files/protein_DNA_affinity.csv', skiprows=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Protein ',\n",
       " 'Specific Binding Constants',\n",
       " 'Non-Specific Binding Constants ',\n",
       " 'Temperature',\n",
       " 'Specific Free energy (presuming room temperature)',\n",
       " 'Specific Free energy (temperature corrected value)',\n",
       " 'Nonspecific Free Energy (presuming room temperature)',\n",
       " 'Nonspecific Free Energy (temperature corrected value)',\n",
       " 'Unnamed: 8']"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# List of column names \n",
    "list(protein_DNA_dfs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Protein</th>\n",
       "      <th>Specific Binding Constants</th>\n",
       "      <th>Non-Specific Binding Constants</th>\n",
       "      <th>Temperature</th>\n",
       "      <th>Specific Free energy (presuming room temperature)</th>\n",
       "      <th>Specific Free energy (temperature corrected value)</th>\n",
       "      <th>Nonspecific Free Energy (presuming room temperature)</th>\n",
       "      <th>Nonspecific Free Energy (temperature corrected value)</th>\n",
       "      <th>Unnamed: 8</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>λ Cro</td>\n",
       "      <td>5.000000e+11</td>\n",
       "      <td>680000.0</td>\n",
       "      <td>273.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>14.62</td>\n",
       "      <td>7.96</td>\n",
       "      <td>7.29</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-</td>\n",
       "      <td>1.200000e+11</td>\n",
       "      <td>NaN</td>\n",
       "      <td>273.0</td>\n",
       "      <td>15.1</td>\n",
       "      <td>13.85</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-</td>\n",
       "      <td>8.300000e+09</td>\n",
       "      <td>NaN</td>\n",
       "      <td>273.0</td>\n",
       "      <td>13.5</td>\n",
       "      <td>12.40</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>λ CcI repressor</td>\n",
       "      <td>1.300000e+09</td>\n",
       "      <td>12000.0</td>\n",
       "      <td>295.0</td>\n",
       "      <td>12.4</td>\n",
       "      <td>12.31</td>\n",
       "      <td>5.57</td>\n",
       "      <td>5.51</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-</td>\n",
       "      <td>8.300000e+09</td>\n",
       "      <td>NaN</td>\n",
       "      <td>293.0</td>\n",
       "      <td>13.5</td>\n",
       "      <td>13.31</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          Protein   Specific Binding Constants  \\\n",
       "0            λ Cro                5.000000e+11   \n",
       "1               -                 1.200000e+11   \n",
       "2               -                 8.300000e+09   \n",
       "3  λ CcI repressor                1.300000e+09   \n",
       "4               -                 8.300000e+09   \n",
       "\n",
       "   Non-Specific Binding Constants   Temperature  \\\n",
       "0                         680000.0        273.0   \n",
       "1                              NaN        273.0   \n",
       "2                              NaN        273.0   \n",
       "3                          12000.0        295.0   \n",
       "4                              NaN        293.0   \n",
       "\n",
       "   Specific Free energy (presuming room temperature)  \\\n",
       "0                                               16.0   \n",
       "1                                               15.1   \n",
       "2                                               13.5   \n",
       "3                                               12.4   \n",
       "4                                               13.5   \n",
       "\n",
       "   Specific Free energy (temperature corrected value)  \\\n",
       "0                                              14.62    \n",
       "1                                              13.85    \n",
       "2                                              12.40    \n",
       "3                                              12.31    \n",
       "4                                              13.31    \n",
       "\n",
       "   Nonspecific Free Energy (presuming room temperature)  \\\n",
       "0                                               7.96      \n",
       "1                                                NaN      \n",
       "2                                                NaN      \n",
       "3                                               5.57      \n",
       "4                                                NaN      \n",
       "\n",
       "   Nonspecific Free Energy (temperature corrected value)  Unnamed: 8  \n",
       "0                                               7.29             NaN  \n",
       "1                                                NaN             NaN  \n",
       "2                                                NaN             NaN  \n",
       "3                                               5.51             NaN  \n",
       "4                                                NaN             NaN  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sample Data \n",
    "protein_DNA_dfs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Protein</th>\n",
       "      <th>Specific Binding Constants</th>\n",
       "      <th>Non-Specific Binding Constants</th>\n",
       "      <th>Temperature</th>\n",
       "      <th>Specific Free energy (presuming room temperature)</th>\n",
       "      <th>Specific Free energy (temperature corrected value)</th>\n",
       "      <th>Nonspecific Free Energy (presuming room temperature)</th>\n",
       "      <th>Nonspecific Free Energy (temperature corrected value)</th>\n",
       "      <th>Unnamed: 8</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>λ Cro</td>\n",
       "      <td>5.000000e+11</td>\n",
       "      <td>680000.0</td>\n",
       "      <td>273.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>14.62</td>\n",
       "      <td>7.96</td>\n",
       "      <td>7.29</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-</td>\n",
       "      <td>1.200000e+11</td>\n",
       "      <td>NaN</td>\n",
       "      <td>273.0</td>\n",
       "      <td>15.1</td>\n",
       "      <td>13.85</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-</td>\n",
       "      <td>8.300000e+09</td>\n",
       "      <td>NaN</td>\n",
       "      <td>273.0</td>\n",
       "      <td>13.5</td>\n",
       "      <td>12.40</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>λ CcI repressor</td>\n",
       "      <td>1.300000e+09</td>\n",
       "      <td>12000.0</td>\n",
       "      <td>295.0</td>\n",
       "      <td>12.4</td>\n",
       "      <td>12.31</td>\n",
       "      <td>5.57</td>\n",
       "      <td>5.51</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-</td>\n",
       "      <td>8.300000e+09</td>\n",
       "      <td>NaN</td>\n",
       "      <td>293.0</td>\n",
       "      <td>13.5</td>\n",
       "      <td>13.31</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          Protein   Specific Binding Constants  \\\n",
       "0            λ Cro                5.000000e+11   \n",
       "1               -                 1.200000e+11   \n",
       "2               -                 8.300000e+09   \n",
       "3  λ CcI repressor                1.300000e+09   \n",
       "4               -                 8.300000e+09   \n",
       "\n",
       "   Non-Specific Binding Constants   Temperature  \\\n",
       "0                         680000.0        273.0   \n",
       "1                              NaN        273.0   \n",
       "2                              NaN        273.0   \n",
       "3                          12000.0        295.0   \n",
       "4                              NaN        293.0   \n",
       "\n",
       "   Specific Free energy (presuming room temperature)  \\\n",
       "0                                               16.0   \n",
       "1                                               15.1   \n",
       "2                                               13.5   \n",
       "3                                               12.4   \n",
       "4                                               13.5   \n",
       "\n",
       "   Specific Free energy (temperature corrected value)  \\\n",
       "0                                              14.62    \n",
       "1                                              13.85    \n",
       "2                                              12.40    \n",
       "3                                              12.31    \n",
       "4                                              13.31    \n",
       "\n",
       "   Nonspecific Free Energy (presuming room temperature)  \\\n",
       "0                                               7.96      \n",
       "1                                                NaN      \n",
       "2                                                NaN      \n",
       "3                                               5.57      \n",
       "4                                                NaN      \n",
       "\n",
       "   Nonspecific Free Energy (temperature corrected value)  Unnamed: 8  \n",
       "0                                               7.29             NaN  \n",
       "1                                                NaN             NaN  \n",
       "2                                                NaN             NaN  \n",
       "3                                               5.51             NaN  \n",
       "4                                                NaN             NaN  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sample Data \n",
    "protein_DNA_dfs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12.916086956521738"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Data table; Mean of the specific free energy in kcal/mol\n",
    "protein_DNA_dfs['Specific Free energy (temperature corrected value)'].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.3808333333333325"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Data table; Mean of the non-specific free energy in kcal/mol\n",
    "protein_DNA_dfs['Nonspecific Free Energy (temperature corrected value)'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (b). Calculating $(N)_{\\text{max}}$, $(NV)_{\\text{max}}$, and Type I inequality\n",
    "\n",
    "We calculate $(N)_{\\text{max}}$, $(NV)_{\\text{max}}$, and the right-hand-side of the \"Type I inequality\" for the TF-DNA system. In contrast to the ssDNA-ssDNA and protein-protein interaction calculations, we use the results for a gendered system. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Relevant Equalities/Inequalities **\n",
    "\n",
    "\\begin{equation}\n",
    "(NV)_{\\text{max}} = \\lambda_{\\mu}^3 \\,e^{\\beta(\\Delta+E_0)}\\qquad \\text{[$T= T_{\\text{I}}'$]},\n",
    "\\end{equation}\n",
    "where $\\lambda_{\\mu} = h/\\sqrt{2\\pi \\mu k_BT}$\n",
    "\n",
    "\\begin{equation}\n",
    "(N)_{\\text{max}} = e^{\\beta \\Delta}\\qquad \\text{[$T= T_{\\text{II}}'$]},\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "N < \\exp\\left[\\frac{3\\Delta}{2E_0} W_{0}\\left(\\frac{2E_0}{3E_{\\mu, V}}\\right)\\right],  \\quad \\text{[Type I inequality, $T_{\\text{I}}< T_{\\text{II}}'$]}\n",
    "\\end{equation}\n",
    "\n",
    "where $E_{\\mu, V} = h^2/2\\pi \\mu V^{2/3}$.\n",
    "\n",
    "**Change in Translational Entropy upon Binding**\n",
    "\n",
    "\\begin{equation}\n",
    " T_{\\text{st}}\\left(S_{\\text{TF-DNA}} -\\,S_{\\text{TF}}-\\,S_{\\text{DNA}}\\right)/N_A =  -T_{\\text{st}} S_{\\text{TF}}/N_A =  - RT_{\\text{st}} \\ln \\left[e^{5/2}\\frac{V_{\\text{st}}}{N_A\\lambda_{\\text{TF}}^3} \\right] = - 16.4 \\text{ kcal/mol},\n",
    "\\end{equation}\n",
    "\n",
    "where we took $S_{\\text{TF-DNA}} = S_{\\text{DNA}}$ since the DNA (binding site) is fixed compared to the motile TF. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-16.355974454934454"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Computing translational entropy contribution to free energy\n",
    "\n",
    "# We use SI units throughout, except for energies which are defined in kcal/mol\n",
    "\n",
    "Tphys = 310.15 # physiological temperature\n",
    "RTphys = 0.001987*Tphys # Gas constant times temperature\n",
    "V_stan = 1e-3 # 1 liter in meters cubed\n",
    "kB = 1.38*10**(-23) # Boltzmann's constant\n",
    "mass = 1.062745*10**(-22) # 64 kDa; mass of protein\n",
    "Na = 6.022*10**23 # Avogadro's number\n",
    "h = 6.626*10**(-34) # Planck's constant\n",
    "lam_0 = h/np.sqrt(2*np.pi*mass*kB*Tphys) # thermal deBroglie wavelength\n",
    "\n",
    "-RTphys*np.log(np.exp(5/2)*V_stan/(Na*lam_0**3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Parameter Estimation**\n",
    "\n",
    "From the specific free energy of binding we find:\n",
    "\n",
    "$$ \\delta E = \\delta G + T \\delta S = -12.9 \\text{ kcal/mol} - 16.4 \\text{ kcal/mol} = -29.3  \\text{ kcal/mol} = -(E_0+\\Delta),$$\n",
    "\n",
    "and from the non-specific free energy of binding we have $$-\\Delta = -6.4\\text{ kcal/mol}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(N)max = 3.24E+04\n",
      "(NV)max= 2.68E+01\n",
      "exp((3Δ/2E0) W0(2E0/3EV))= 2.19E+05\n"
     ]
    }
   ],
   "source": [
    "## Calculation (N)max, (NV)max, and Type I inequality\n",
    "\n",
    "# System paramaters for DNA system\n",
    "kB = 1.38*10**(-23) # Boltzmann's constant\n",
    "mass = 1.062745*10**(-22) # 64 kDa; mass of protein\n",
    "h = 6.626*10**(-34) # Planck's constant\n",
    "Na = 6.022*10**23 # Avogadro's number\n",
    "\n",
    "E0 = 29.3-6.4\n",
    "Del = 6.4\n",
    "Tphys = 310.15 # physiological temperature\n",
    "RTphys = 0.001987*Tphys # Gas constant times temperature\n",
    "Vol = 1e-18 # we assume a volume of (1 micrometer)^3\n",
    "\n",
    "lam_0_phys = h/np.sqrt(2*np.pi*mass*kB*Tphys) # thermal deBroglie wavelength at phys temperature\n",
    "\n",
    "Ev = h**2/(2*np.pi*mass*Vol**(2/3))*(2.388e-4)*Na # Ev written in kcal/mol\n",
    "\n",
    "Nmax = np.exp(Del/RTphys)\n",
    "NVmax = lam_0_phys**3*np.exp((E0+Del)/RTphys)/(1e-18)\n",
    "Type1_ineq = np.exp((3*Del)/(2*E0)*lambertw(2*E0/(3*Ev))).real \n",
    "\n",
    "print(\"(N)max =\", '%.2E' % Decimal(Nmax))\n",
    "print(\"(NV)max=\", '%.2E' % Decimal(NVmax))\n",
    "print(\"exp((3Δ/2E0) W0(2E0/3EV))=\", '%.2E' % Decimal(Type1_ineq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----\n",
    "\n",
    "------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### III. Protein-protein interaction calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (a). Free Energy Estimate\n",
    "\n",
    "\n",
    "We estimate the specific free energy of binding for protein-protein intereactions for some representative interactions\n",
    "\n",
    "Data is taken from the \"Affinity Benchmark 2.0 Values\" here https://bmm.crick.ac.uk/~bmmadmin/Affinity/\n",
    "\n",
    "* Kastritis, P. L., Moal, I. H., Hwang, H., Weng, Z., Bates, P. A., Bonvin, A. M., & Janin, J. (2011). A structure‐based benchmark for protein–protein binding affinity. *Protein Science*, 20(3), 482-491.\n",
    "\n",
    "Difference between non-specific an specific binding (i.e., $(\\delta G)_{\\text{non-spec.}}- (\\delta G)_{\\text{spec.}}$) is found to be $12.5\\,RT$ (last paragraph of pg. 5) or $7.7$ kcal/mol for $T = 310.15$ K. \n",
    "\n",
    "* Zhang, J., Maslov, S., & Shakhnovich, E. I. (2008). Constraints imposed by non‐functional protein–protein interactions on gene expression and proteome size. *Molecular systems biology*, 4(1), 210."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 12 TFs and their specific and non-specific binding energies to DNA\n",
    "'''\n",
    "Data from:\n",
    "\n",
    "Protein-Protein Affinity Database: https://bmm.crick.ac.uk/~bmmadmin/Affinity/\n",
    "\n",
    "associated with paper \n",
    "\n",
    "Kastritis, Panagiotis L., et al. \"A structure‐based benchmark \n",
    "for protein–protein binding affinity.\" Protein Science 20.3 (2011): 482-491.\n",
    "'''\n",
    "protein_protein_dfs = pd.read_csv('data_files/protein_protein_affinity.csv', skiprows=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Complex PDB',\n",
       " 'Type',\n",
       " 'Unbound PDB',\n",
       " 'Protein A',\n",
       " 'Unbound PDB.1',\n",
       " 'Protein B',\n",
       " 'Pubmed',\n",
       " 'Kd (M)',\n",
       " 'dG',\n",
       " 'I-RMSD',\n",
       " 'DASA ',\n",
       " 'Method',\n",
       " 'Final Comments',\n",
       " 'Temperature (C)',\n",
       " 'pH',\n",
       " 'Buffer/Ionic Strength',\n",
       " 'Cofactors (bound)',\n",
       " 'Cofactors (unbound)',\n",
       " 'Cofactors (binding)',\n",
       " 'Extra Data; NB) kon in mol^-1 s^-1, koff in s^-1, dH in kcal mol^-1 dS and dCp in cal mol^-1 K^-1',\n",
       " 'Corroborating Data',\n",
       " 'Exact finding in: ']"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# List of column names \n",
    "list(protein_protein_dfs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Complex PDB</th>\n",
       "      <th>Type</th>\n",
       "      <th>Unbound PDB</th>\n",
       "      <th>Protein A</th>\n",
       "      <th>Unbound PDB.1</th>\n",
       "      <th>Protein B</th>\n",
       "      <th>Pubmed</th>\n",
       "      <th>Kd (M)</th>\n",
       "      <th>dG</th>\n",
       "      <th>I-RMSD</th>\n",
       "      <th>...</th>\n",
       "      <th>Final Comments</th>\n",
       "      <th>Temperature (C)</th>\n",
       "      <th>pH</th>\n",
       "      <th>Buffer/Ionic Strength</th>\n",
       "      <th>Cofactors (bound)</th>\n",
       "      <th>Cofactors (unbound)</th>\n",
       "      <th>Cofactors (binding)</th>\n",
       "      <th>Extra Data; NB) kon in mol^-1 s^-1, koff in s^-1, dH in kcal mol^-1 dS and dCp in cal mol^-1 K^-1</th>\n",
       "      <th>Corroborating Data</th>\n",
       "      <th>Exact finding in:</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1A2K_C:AB</td>\n",
       "      <td>OG</td>\n",
       "      <td>1QG4_A</td>\n",
       "      <td>Ran GTPase-GDP</td>\n",
       "      <td>1OUN_AB</td>\n",
       "      <td>Nuclear transport factor 2</td>\n",
       "      <td>10681579.0</td>\n",
       "      <td>0.00000015</td>\n",
       "      <td>-9.31</td>\n",
       "      <td>1.11</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>25</td>\n",
       "      <td>7.5</td>\n",
       "      <td>1</td>\n",
       "      <td>GDP/Mg</td>\n",
       "      <td>GDP/Mg</td>\n",
       "      <td>GDP/Mg</td>\n",
       "      <td>Microtitre plates: Kd = 2.4e-7.  Fluorescence:...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Fig. 1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1ACB_E:I</td>\n",
       "      <td>EI</td>\n",
       "      <td>4CHA_ABC</td>\n",
       "      <td>Chymotrypsin</td>\n",
       "      <td>1EGL_A</td>\n",
       "      <td>Eglin C</td>\n",
       "      <td>3071573.0</td>\n",
       "      <td>0.0000000002</td>\n",
       "      <td>-13.05</td>\n",
       "      <td>1.08</td>\n",
       "      <td>...</td>\n",
       "      <td>calculated from Ka</td>\n",
       "      <td>21</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>DG(STP) = -13.0, dH(STP) = +2.0, dS(STP) = +51...</td>\n",
       "      <td>PMID 1741752:  Kd = 4e-13, (50mM Hepes, 100mM ...</td>\n",
       "      <td>Table 1, calculated from Ka; PMID 1741752: p.31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1AHW_AB:C</td>\n",
       "      <td>A</td>\n",
       "      <td>1FGN_LH</td>\n",
       "      <td>Fab 5g9</td>\n",
       "      <td>1TFH_A</td>\n",
       "      <td>Tissue factor</td>\n",
       "      <td>9480775.0</td>\n",
       "      <td>0.0000000034</td>\n",
       "      <td>-11.55</td>\n",
       "      <td>0.69</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>ambient</td>\n",
       "      <td>not stated</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Table 4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1AK4_A:D</td>\n",
       "      <td>OX</td>\n",
       "      <td>2CPL_A</td>\n",
       "      <td>Cyclophilin</td>\n",
       "      <td>4J93_A</td>\n",
       "      <td>HIV capsid</td>\n",
       "      <td>9223641.0</td>\n",
       "      <td>0.000016</td>\n",
       "      <td>-6.43</td>\n",
       "      <td>1.33</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>20</td>\n",
       "      <td>6.5</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SPR (same conditions): Kd = 1.8e-5.  Equillibr...</td>\n",
       "      <td>PMID 19767750 (ITC):  Kd = 5e-5, dCp = -285 dH...</td>\n",
       "      <td>Fig. 4; PMID 19767750: Fig.2, PMID 8980234: p....</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1AKJ_AB:DE</td>\n",
       "      <td>OX</td>\n",
       "      <td>2CLR_DE</td>\n",
       "      <td>MHC class 1 HLA-A2</td>\n",
       "      <td>1CD8_AB</td>\n",
       "      <td>T-cell CD8 coreceptor</td>\n",
       "      <td>10072074.0</td>\n",
       "      <td>0.000126</td>\n",
       "      <td>-5.32</td>\n",
       "      <td>1.14</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>25</td>\n",
       "      <td>7.4</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Kon ? 140000, koff ? 18 s?1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Table 1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 22 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  Complex PDB Type Unbound PDB           Protein A Unbound PDB.1  \\\n",
       "0   1A2K_C:AB   OG      1QG4_A      Ran GTPase-GDP       1OUN_AB   \n",
       "1    1ACB_E:I   EI    4CHA_ABC        Chymotrypsin        1EGL_A   \n",
       "2   1AHW_AB:C    A     1FGN_LH             Fab 5g9        1TFH_A   \n",
       "3    1AK4_A:D   OX      2CPL_A         Cyclophilin        4J93_A   \n",
       "4  1AKJ_AB:DE   OX     2CLR_DE  MHC class 1 HLA-A2       1CD8_AB   \n",
       "\n",
       "                    Protein B      Pubmed        Kd (M)     dG  I-RMSD  ...  \\\n",
       "0  Nuclear transport factor 2  10681579.0    0.00000015  -9.31    1.11  ...   \n",
       "1                     Eglin C   3071573.0  0.0000000002 -13.05    1.08  ...   \n",
       "2               Tissue factor   9480775.0  0.0000000034 -11.55    0.69  ...   \n",
       "3                  HIV capsid   9223641.0      0.000016  -6.43    1.33  ...   \n",
       "4       T-cell CD8 coreceptor  10072074.0      0.000126  -5.32    1.14  ...   \n",
       "\n",
       "       Final Comments Temperature (C)          pH Buffer/Ionic Strength  \\\n",
       "0                 NaN              25         7.5                     1   \n",
       "1  calculated from Ka              21           8                     1   \n",
       "2                 NaN         ambient  not stated                     1   \n",
       "3                 NaN              20         6.5                     1   \n",
       "4                 NaN              25         7.4                     1   \n",
       "\n",
       "  Cofactors (bound)  Cofactors (unbound) Cofactors (binding)  \\\n",
       "0            GDP/Mg               GDP/Mg              GDP/Mg   \n",
       "1               NaN                  NaN                 NaN   \n",
       "2               NaN                  NaN                 NaN   \n",
       "3               NaN                  NaN                 NaN   \n",
       "4               NaN                  NaN                 NaN   \n",
       "\n",
       "  Extra Data; NB) kon in mol^-1 s^-1, koff in s^-1, dH in kcal mol^-1 dS and dCp in cal mol^-1 K^-1  \\\n",
       "0  Microtitre plates: Kd = 2.4e-7.  Fluorescence:...                                                  \n",
       "1  DG(STP) = -13.0, dH(STP) = +2.0, dS(STP) = +51...                                                  \n",
       "2                                                NaN                                                  \n",
       "3  SPR (same conditions): Kd = 1.8e-5.  Equillibr...                                                  \n",
       "4                        Kon ? 140000, koff ? 18 s?1                                                  \n",
       "\n",
       "                                  Corroborating Data  \\\n",
       "0                                                NaN   \n",
       "1  PMID 1741752:  Kd = 4e-13, (50mM Hepes, 100mM ...   \n",
       "2                                                NaN   \n",
       "3  PMID 19767750 (ITC):  Kd = 5e-5, dCp = -285 dH...   \n",
       "4                                                NaN   \n",
       "\n",
       "                                  Exact finding in:   \n",
       "0                                             Fig. 1  \n",
       "1    Table 1, calculated from Ka; PMID 1741752: p.31  \n",
       "2                                            Table 4  \n",
       "3  Fig. 4; PMID 19767750: Fig.2, PMID 8980234: p....  \n",
       "4                                            Table 1  \n",
       "\n",
       "[5 rows x 22 columns]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sample Data \n",
    "protein_protein_dfs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-10.877608695652173"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Data table; Mean of the binding free energy in kcal/mol\n",
    "protein_protein_dfs['dG'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (b). Calculating $(N)_{\\text{max}}$, $(NV)_{\\text{max}}$, and Type I inequality\n",
    "\n",
    "We calculate $(N)_{\\text{max}}$, $(NV)_{\\text{max}}$, and the right-hand-side of the \"Type I inequality\" for protein-protein systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Relevant Inequalities **\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "(NV)_{\\text{max}} = \\sqrt{2} \\lambda_{0}^3 e^{\\beta(\\Delta+E_0)}, \\qquad \\text{[$T= T_{\\text{I}}$]}\n",
    "\\end{equation}\n",
    "where $\\lambda_{0} = h/\\sqrt{2\\pi m_0 k_BT}$\n",
    "\n",
    "\\begin{equation}\n",
    "(N)_{\\text{max}} = e^{\\beta \\Delta}/2, \\qquad \\text{[$T= T_{\\text{II}}$]}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "2N < \\exp\\left[\\frac{3\\Delta}{2E_0} W_{0}\\left(\\frac{E_0}{3E_V}\\right)\\right], \\qquad \\text{[Type I inequality, $T_{\\text{I}}< T_{\\text{II}}$]}\n",
    "\\end{equation}\n",
    "\n",
    "where $E_{ V} = h^2/2\\pi m_0 V^{2/3}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Change in Translational Entropy upon Binding**\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    " T_{\\text{st}}\\left(S_{\\text{protein-protein}} - 2 S_{\\text{protein}}\\right)/N_A = - R T_{\\text{st}} \\ln \\left[e^{5/2}\\frac{V_{\\text{st}}}{2\\sqrt{2}\\,N_A\\lambda_{0}^3} \\right] = - 15.7 \\text{ kcal/mol}.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-15.71522776294444"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Computing translational entropy contribution to free energy\n",
    "\n",
    "# We use SI units throughout, except for energies which are defined in kcal/mol\n",
    "\n",
    "Tphys = 310.15 # physiological temperature\n",
    "RTphys = 0.001987*Tphys # Gas constant times temperature\n",
    "V_stan = 1e-3 # 1 liter in meters cubed\n",
    "kB = 1.38*10**(-23) # Boltzmann's constant\n",
    "mass = 1.062745*10**(-22) # 64 kDa; mass of protein\n",
    "Na = 6.022*10**23 # Avogadro's number\n",
    "h = 6.626*10**(-34) # Planck's constant\n",
    "lam_0 = h/np.sqrt(2*np.pi*mass*kB*Tphys) # thermal deBroglie wavelength\n",
    "\n",
    "-RTphys*np.log(np.exp(5/2)*V_stan/(2*np.sqrt(2)*Na*lam_0**3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Parameter Estimation**\n",
    "\n",
    "From the free energy of binding (presumed to be specific) we find:\n",
    "\n",
    "$$ \\delta E = \\delta G + T \\delta S = -10.9 \\text{ kcal/mol} - 15.7 \\text{ kcal/mol} = -26.6  \\text{ kcal/mol} = -(E_0+\\Delta)$$\n",
    "\n",
    "and from Zhang *et.al.* refernce we have $$-\\Delta = -7.7\\text{ kcal/mol}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(N)max = 1.33E+05\n",
      "(NV)max= 4.74E-01\n",
      "exp((3Δ/2E0) W0(E0/3EV))= 3.62E+07\n"
     ]
    }
   ],
   "source": [
    "## Calculation (N)max, (NV)max, and Type I inequality\n",
    "\n",
    "# System paramaters for DNA system\n",
    "kB = 1.38*10**(-23) # Boltzmann's constant\n",
    "mass = 1.062745*10**(-22) # 64 kDa; mass of protein\n",
    "h = 6.626*10**(-34) # Planck's constant\n",
    "Na = 6.022*10**23 # Avogadro's number\n",
    "\n",
    "E0 = 26.6-7.7\n",
    "Del = 7.7\n",
    "Tphys = 310.15 # physiological temperature\n",
    "RTphys = 0.001987*Tphys # Gas constant times temperature\n",
    "Vol = 1e-18 # we assume a volume of (1 micrometer)^3\n",
    "\n",
    "lam_0_phys = h/np.sqrt(2*np.pi*mass*kB*Tphys) # thermal deBroglie wavelength at phys temperature\n",
    "\n",
    "Ev = h**2/(2*np.pi*mass*Vol**(2/3))*(2.388e-4)*Na # Ev written in kcal/mol\n",
    "\n",
    "Nmax = np.exp(Del/RTphys)/2\n",
    "NVmax = np.sqrt(2)*lam_0_phys**3*np.exp((E0+Del)/RTphys)/(1e-18)\n",
    "Type1_ineq = np.exp((3*Del)/(2*E0)*lambertw(E0/(3*Ev))).real \n",
    "\n",
    "\n",
    "print(\"(N)max =\", '%.2E' % Decimal(Nmax))\n",
    "print(\"(NV)max=\", '%.2E' % Decimal(NVmax))\n",
    "print(\"exp((3Δ/2E0) W0(E0/3EV))=\", '%.2E' % Decimal(Type1_ineq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed Time: 3.276 sec\n"
     ]
    }
   ],
   "source": [
    "print('Elapsed Time: %.3f sec' % (time.time()-start_clock))"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}