index.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\var}[3]{{{#1}_\\text{#2}^\\text{#3}}}$\n",
    "\n",
    "$\\newcommand{\\pa}{p_\\text{A}}$\n",
    "\n",
    "$\\newcommand{\\pb}{p_\\text{B}}$\n",
    "\n",
    "$\\newcommand{\\wa}{ω_\\text{A}}$\n",
    "\n",
    "$\\newcommand{\\wb}{ω_\\text{B}}$\n",
    "\n",
    "$\\newcommand{\\kab}{k_\\text{AB}}$ $\\newcommand{\\kba}{k_\\text{BA}}$\n",
    "\n",
    "$\\newcommand{\\kex}{k_\\text{ex}}$ $\\newcommand{\\ddt}{\\frac{d}{dt}}$\n",
    "\n",
    "$\\newcommand{\\kon}{k_\\text{on}}$ $\\newcommand{\\koff}{k_\\text{off}}$\n",
    "\n",
    "$\\newcommand{\\nuf}{\\nu_\\text{free}}$\n",
    "\n",
    "$\\newcommand{\\nub}{\\nu_\\text{bound}}$\n",
    "\n",
    "$\\newcommand{\\num}{\\nu_\\text{measure}}$\n",
    "\n",
    "$\\newcommand{\\Kd}{K_\\text{d}}$\n",
    "\n",
    "---\n",
    "\n",
    "# Chemical Exchange in NMR Spectroscopy\n",
    "\n",
    "**Guillaume Bouvignies**, **Bruno Kieffer**, **Catherine Bougault**, **Nathalie Sibille**\n",
    "\n",
    "_RéNaFoBiS, Oleron, 2024_\n",
    "\n",
    "---\n",
    "\n",
    "## 1. The Model\n",
    "\n",
    "In this section, we examine a protein undergoing chemical exchange between two states, **A** and **B**, each characterized by distinct chemical environments. The exchange process is represented by the equation:\n",
    "\n",
    "$$ A \\underset{\\kba}{\\stackrel{\\kab}{\\rightleftharpoons}} B \\tag{1} $$\n",
    "\n",
    "Here, $\\kab$ and $\\kba$ are the forward and reverse rate constants, respectively, both of which are first-order. The sum of these rate constants, $\\kex = \\kab + \\kba$, represents the overall exchange rate. The fractional populations of states **A** and **B** are denoted as $\\pa$ and $\\pb$, respectively, with the constraint that $\\pa + \\pb = 1$. The lifetimes of states **A** and **B** are defined as $\\tau_A = 1/\\kab$ and $\\tau_B = 1/\\kba$. The timescale of the exchange process is defined as $\\tau_{ex} = 1 / \\kex$.\n",
    "\n",
    "In the subsequent discussion, we focus on a specific nucleus within the protein. For this nucleus, the chemical shift (measured in rad·s<sup>-1</sup>) is $\\wa$ in state **A** and $\\wb$ in state **B**, with a difference of $\\Delta \\omega_{AB} = \\wb - \\wa$.\n",
    "\n",
    "## 2. The Kinetic Equations\n",
    "\n",
    "The differential equation system for the first-order two-site exchange reaction, as given in Equation (1), can be formulated using the chemical reaction rate law:\n",
    "\n",
    "$$ \\ddt[A] = - \\kab[A] + \\kba[B] \\tag{2} $$\n",
    "\n",
    "$$ \\ddt[B] = - \\kba[B] + \\kab[A] \\tag{3} $$\n",
    "\n",
    "These equations can be expressed in matrix form as:\n",
    "\n",
    "$$\n",
    "\\ddt \\begin{pmatrix} [A] \\\\ [B] \\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "    -\\kab & \\kba \\\\\n",
    "    \\kab & -\\kba \\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix} [A] \\\\ [B] \\end{pmatrix}\n",
    "\\tag{4}\n",
    "$$\n",
    "\n",
    "At equilibrium, we have:\n",
    "\n",
    "$$\n",
    "\\ddt[A]=\\ddt[B]=0 \\implies \\frac{\\var{[A]}{eq}{}}{\\var{[B]}{eq}{}}=\\frac{\\kba}{\\kab} \\tag{5}\n",
    "$$\n",
    "\n",
    "which leads to:\n",
    "\n",
    "$$\n",
    "\\pb = \\frac{\\var{[B]}{eq}{}}{\\var{[A]}{eq}{}+\\var{[B]}{eq}{}} = \\frac{\\kab}{\\kab+\\kba} \\tag{6}\n",
    "$$\n",
    "\n",
    "## 3. The Bloch Equations\n",
    "\n",
    "The Bloch equations describe the evolution of magnetization (in the rotating frame) with the $\\text{B}_0$ field aligned in the z-direction. These can be represented in matrix form as:\n",
    "\n",
    "$$\n",
    "\\ddt \\begin{pmatrix} M_\\text{x} \\\\ M_\\text{y} \\\\ M_\\text{z} \\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "    -R_2 & -Ω & ω_1\\sin\\phi \\\\\n",
    "    Ω & -R_2 & -ω_1\\cos\\phi \\\\\n",
    "    -ω_1\\sin\\phi & ω_1\\cos\\phi & -R_1 \\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix} M_\\text{x} \\\\ M_\\text{y} \\\\ M_\\text{z} \\end{pmatrix} +\n",
    "\\begin{pmatrix} 0 \\\\ 0 \\\\ R_1 \\var{M}{eq}{} \\end{pmatrix}\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "Here, $R_1$ and $R_2$ represent the rates at which the longitudinal and transverse magnetization return to equilibrium (in s<sup>-1</sup>), $Ω$ is the offset resonance frequency (in rad·s<sup>-1</sup>), $ω_1$ and $\\phi$ are the strength (in rad·s<sup>-1</sup>) and phase of the applied RF field, respectively. $M_{i}$ denotes the magnetization component along the $i \\in \\{x,y,z\\}$ axis, and $\\var{M}{eq}{}$ refers to the equilibrium magnetization.\n",
    "\n",
    "To express the Bloch equations in a homogeneous form, we augment the magnetization component vector with the identity operator $E$:\n",
    "\n",
    "$$\n",
    "\\ddt\n",
    "\\begin{pmatrix} E \\\\ M_\\text{x} \\\\ M_\\text{y} \\\\ M_\\text{z} \\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "    0 & 0 & 0 & 0 \\\\\n",
    "    0 & -R_2 & -Ω & ω_1\\sin\\phi \\\\\n",
    "    0 & Ω & -R_2 & -ω_1\\cos\\phi \\\\\n",
    "    R_1\\var{M}{eq}{} & -ω_1\\sin\\phi & ω_1\\cos\\phi & -R_1 \\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix} E \\\\ M_\\text{x} \\\\ M_\\text{y} \\\\ M_\\text{z} \\end{pmatrix}\n",
    "\\tag{8}\n",
    "$$\n",
    "\n",
    "## 4. The Bloch-McConnell Equations\n",
    "\n",
    "In a seminal 1958 article, McConnell observed the analogy between magnetization and concentration, demonstrating that the Bloch equations can be extended to include the effects of chemical exchange. The resulting generalized equation can be expressed in matrix form as follows:\n",
    "\n",
    "$$\n",
    "\\ddt \\begin{pmatrix} E \\\\ M_\\text{x}^\\text{A} \\\\ M_\\text{y}^\\text{A} \\\\ M_\\text{z}^\\text{A} \\\\ M_\\text{x}^\\text{B} \\\\ M_\\text{y}^\\text{B} \\\\ M_\\text{z}^\\text{B} \\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "    0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n",
    "    0 & -R_{2}^\\text{A} - \\kab & -Ω^\\text{A} & ω_1\\sin\\phi & \\kba & 0 & 0 \\\\\n",
    "    0 & Ω^\\text{A} & -R_{2}^\\text{A} - \\kab & -ω_1\\cos\\phi & 0 & \\kba & 0 \\\\\n",
    "    R_{1}^\\text{A}M_{eq}^\\text{A} & -ω_1\\sin\\phi & ω_1\\cos\\phi & -R_{1}^\\text{A} - \\kab & 0 & 0 & \\kba \\\\\n",
    "    0 & \\kab & 0 & 0 & -R_{2}^\\text{B} - \\kba & -Ω^\\text{B} & ω_1\\sin\\phi \\\\\n",
    "    0 & 0 & \\kab & 0 & -Ω^\\text{B} & -R_{2}^\\text{B} - \\kba & -ω_1\\cos\\phi \\\\\n",
    "    R_{1}^\\text{B}M_{eq}^\\text{B} & 0 & 0 & \\kab & -ω_1\\sin\\phi & ω_1\\cos\\phi & -R_{1}^\\text{B} - \\kba \\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix} E \\\\ M_\\text{x}^\\text{A} \\\\ M_\\text{y}^\\text{A} \\\\ M_\\text{z}^\\text{A} \\\\ M_\\text{x}^\\text{B} \\\\ M_\\text{y}^\\text{B} \\\\ M_\\text{z}^\\text{B} \\end{pmatrix} \\tag{9}\n",
    "$$\n",
    "\n",
    "In Equation (9), state-specific parameters are denoted by the superscript $i \\in \\{ A, B \\}$. It should be noted that Equation (9) corresponds to the scenario of a single isolated spin-½ undergoing exchange between two chemical states. However, the same derivation can be applied to more complex spin systems and kinetic models.\n",
    "\n",
    "The solution to the homogeneous first-order differential equation:\n",
    "\n",
    "$$ \\ddt\\mathbf{M}(t) = \\mathbf{L} \\cdot \\mathbf{M}(t) \\tag{10} $$\n",
    "\n",
    "is given by:\n",
    "\n",
    "$$\n",
    "\\mathbf{M}(t+\\Delta t) = \\exp{\\left(\\mathbf{L}\\Delta t \\right)} \\cdot \\mathbf{M}(t) \\tag{11}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{L}$ and $\\mathbf{M}$ correspond to the matrix and vector, respectively, in Equations (5), (8), and (9).\n",
    "\n",
    "To simulate the magnetization trajectory during an NMR experiment in the presence of chemical exchange, the propagator $\\exp{\\left( \\mathbf{L} \\Delta t \\right)}$ needs to be computed for each individual part of the experiment and multiplied successively with the initial magnetization vector.\n",
    "\n",
    "## 5. The 1D NMR Spectrum in Presence of Two-Site Exchange\n",
    "\n",
    "To simulate the spectra obtained after a 90° pulse, we only need to consider the transverse magnetization component $M^+ = M_\\text{x} + iM_\\text{y}$. The Bloch-McConnell equations in the laboratory frame can then be simplified as:\n",
    "\n",
    "$$\n",
    "\\ddt \\begin{pmatrix} M_\\text{A}^+ \\\\ M_\\text{B}^+ \\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "    -R_{2}^\\text{A} + i\\wa - \\kab & \\kba \\\\\n",
    "    \\kab & -R_{2}^\\text{B} + i\\wb - \\kba \\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix} M_\\text{A}^+ \\\\ M_\\text{B}^+ \\end{pmatrix} \\tag{12}\n",
    "$$\n",
    "\n",
    "The solution to the equation is given by:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix} M_\\text{A}^+(t) \\\\ M_\\text{B}^+(t) \\end{pmatrix} =\n",
    "\\begin{pmatrix} a_{11} & a_{12} \\\\ a_{21} & a_{22} \\end{pmatrix}\n",
    "\\begin{pmatrix} M_\\text{A}^+(0) \\\\ M_\\text{B}^+(0) \\end{pmatrix} \\tag{13}\n",
    "$$\n",
    "\n",
    "where the matrix elements are given by:\n",
    "\n",
    "$$\n",
    "a_{11}(t) = \\frac{1}{2}\n",
    "\\left[\n",
    "    \\left( 1 - \\frac\n",
    "        {i\\Delta \\omega - \\Delta R_2 + \\kex (\\pb - \\pa)}\n",
    "        {λ_+ - λ_-}\n",
    "    \\right)\n",
    "    e^{-λ_-t} +\n",
    "    \\left( 1 + \\frac\n",
    "        {i\\Delta \\omega - \\Delta R_2 + \\kex (\\pb - \\pa)}\n",
    "        {λ_+ - λ_-}\n",
    "    \\right)\n",
    "    e^{-λ_+t}\n",
    "\\right] \\tag{14}\n",
    "$$\n",
    "\n",
    "$$\n",
    "a_{22}(t) = \\frac{1}{2}\n",
    "\\left[\n",
    "    \\left( 1 + \\frac\n",
    "        {i\\Delta \\omega - \\Delta R_2 + \\kex (\\pb - \\pa)}\n",
    "        {λ_+ - λ_-}\n",
    "    \\right)\n",
    "    e^{-λ_-t} +\n",
    "    \\left( 1 - \\frac\n",
    "        {i\\Delta \\omega - \\Delta R_2 + \\kex (\\pb - \\pa)}\n",
    "        {λ_+ - λ_-}\n",
    "    \\right)\n",
    "    e^{-λ_+t}\n",
    "\\right] \\tag{15}\n",
    "$$\n",
    "\n",
    "$$\n",
    "a_{12}(t) = \\frac{\\kex \\pa}{λ_+ - λ_-} \\left( e^{-λ_-t} - e^{-λ_+t} \\right) \\tag{16}\n",
    "$$\n",
    "\n",
    "$$\n",
    "a_{21}(t) = \\frac{\\kex \\pb}{λ_+ - λ_-} \\left( e^{-λ_-t} - e^{-λ_+t} \\right) \\tag{17}\n",
    "$$\n",
    "\n",
    "where $\\Delta R_2 = R_{2}^\\text{B} - R_{2}^\\text{A}$ and the eigenvalues of the Bloch-McConnell equations are given by:\n",
    "\n",
    "$$\n",
    "λ_± = \\frac{1}{2} \\left[\n",
    "    -i (\\wa + \\wb) + R_{2}^\\text{A} + R_{2}^\\text{B} + \\kex \\pm\n",
    "    \\sqrt{\n",
    "        \\left(\n",
    "            i\\Delta \\omega - \\Delta R_2 + \\kex (\\pb - \\pa)\n",
    "        \\right)^2 +\n",
    "        4 \\kex^2 \\pa \\pb\n",
    "    }\n",
    "\\right] \\tag{18}\n",
    "$$\n",
    "\n",
    "## 6. Exchange Regimes\n",
    "\n",
    "<img src=\"Images/chemical_exchange.png\" width=\"60%\">\n",
    "\n",
    "**Figure 1**: Spectra of <sup>15</sup>N nucleus undergoing chemical exchange in an 800 MHz magnet, simulated using the Bloch-McConnell equations. $\\Delta \\delta_{AB} = 3$ ppm and $\\pb = 0.5$ on the left and $\\pb = 0.05$ on the right.\n",
    "\n",
    "**Figure 1** shows a series of spectra calculated using the Bloch-McConnell equations (13) as a function of $\\kex$. The left panel illustrates the scenario where the populations of the two states are equal ($\\pa = \\pb = 0.5$). When the exchange is slow ($\\kex \\ll \\Delta \\omega_{AB}$), two peaks are present. These peaks broaden as the exchange rate increases and eventually merge into a single, very broad peak at $\\kex \\approx \\Delta \\omega_{AB}$, a phenomenon known as coalescence. In the fast exchange regime ($\\kex \\gg \\Delta \\omega_{AB}$), a narrow peak appears at the average position.\n",
    "\n",
    "The more intriguing situation arises when the populations are highly skewed, with $\\pa \\gg \\pb$ (right panel of **Figure 1**). In the slow exchange regime, two peaks are present at positions $ω_{A}$ and $ω_{B}$, with their volumes proportional to the populations. As the exchange rate increases, the peak corresponding to the minor state (**B**) experiences significant line broadening, rapidly becoming invisible. The peak corresponding to the major state (**A**), however, remains visible and broadens as it approaches the coalescence limit, albeit to a lesser extent than resonance **B**. As the exchange becomes faster, the peak narrows and ultimately appears at the position corresponding to the weighted average of the two states.\n",
    "\n",
    "The exchange regime is determined by the ratio of the exchange rate ($\\kex$) to the chemical shift difference between the two states ($\\Delta \\omega_{AB}$), $\\kex / \\Delta \\omega_{AB}$. Thus, one can push the system toward a slower (faster) exchange regime either by decreasing (increasing) $\\kex$, which can be achieved by changing the experimental conditions such as temperature, pH, or ionic strength of the buffer, or by increasing (decreasing) $\\Delta \\omega_{AB}$, which is easily done by running the experiment at a higher (lower) static magnetic field, $\\text{B}_0$.\n",
    "\n",
    "These observations can be further quantitatively described by examining the approximate solutions of the Bloch-McConnell equations for the two limiting cases.\n",
    "\n",
    "### 6.1. Slow Exchange\n",
    "\n",
    "In the slow exchange regime, $\\kex \\ll \\left| i\\Delta \\omega_{AB} + \\Delta R_2 \\right|$, Equation (13) simplifies to:\n",
    "\n",
    "$$\n",
    "M^+(t) = M_\\text{A}^+(0) \\exp\\left[\\left( i\\omega_A - R_2^A - \\kab \\right)t\\right] +\n",
    "M_\\text{B}^+(0) \\exp\\left[\\left( i\\omega_B - R_2^B - \\kba \\right)t\\right] \\tag{19}\n",
    "$$\n",
    "\n",
    "As expected, we observe two resonances at positions corresponding to states **A** and **B**. However, the exchange process adds to the effective relaxation rates of the two peaks, leading to larger linewidths, with $R_{\\text{2,eff}}^\\text{A} = R_2^A + \\kab$ and $R_{\\text{2,eff}}^\\text{B} = R_2^B + \\kba$ for states **A** and **B**, respectively. In the case where $\\pa \\gg \\pb$, i.e., $\\kba \\gg \\kab$, the exchange contribution to the linewidth of peak **B** becomes much larger than that of peak **A**. For example, assuming $\\kex = 1000 \\text{ s}^{-1}$ and $\\pb = 2\\%$, we have $\\kab = 20 \\text{ s}^{-1}$ and $\\kba = 980 \\text{ s}^{-1}$, making peak **B** too broad to be detected in the NMR spectra.\n",
    "\n",
    "### 6.2. Fast Exchange\n",
    "\n",
    "In the fast exchange regime, $\\kex \\gg \\left| i\\Delta \\omega_{AB} + \\Delta R_2 \\right|$, Equation (13) approximates to:\n",
    "\n",
    "$$\n",
    "M^+(t) = M^+(0) \\exp\\left[\\left( i\\overline{\\omega} - \\overline{R_2} -\n",
    "\\frac{\\pa\\pb \\Delta \\omega_{AB}^2}{\\kex} \\right)t\\right] \\tag{20}\n",
    "$$\n",
    "\n",
    "where $\\overline{R_2} = \\pa R_2^A + \\pb R_2^B$ and $\\overline{\\omega} = \\pa \\omega_A + \\pb \\omega_B$ are the population-averaged transverse relaxation rate and chemical shift, respectively. Here, only a single line at the weighted-average position is observed with a linewidth corresponding to the average relaxation rate, plus a contribution due to the exchange process. Interestingly, while the exchange-induced relaxation rate increase is independent of the $\\text{B}_0$ field in the slow exchange regime, it varies quadratically with $\\omega_{AB}$, and thus as $\\text{B}_0^{2}$, in the fast exchange regime.\n",
    "\n",
    "## 7. Simulation of Two-State Exchange\n",
    "\n",
    "Run the following cells to see the effect of different parameters on the NMR line shapes.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "import ipywidgets as widgets\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from simulator import calculate_spectrum, Parameters\n",
    "\n",
    "spectral_width1 = 15.0\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "\n",
    "def update_values(b0, kex, pb, δa, δb, r2a, r2b, yscale):\n",
    "    params = Parameters(b0, pb, kex, δa, δb, r2a, r2b)\n",
    "    coordinates = calculate_spectrum(spectral_width1, params)\n",
    "\n",
    "    ylim = ax.get_ylim()\n",
    "    ax.cla()\n",
    "    ax.plot(coordinates.x, coordinates.y1, \"--\", color=\"C1\")\n",
    "    ax.plot(coordinates.x, coordinates.y2, \"--\", color=\"C2\")\n",
    "    ax.plot(coordinates.x, coordinates.y, color=\"C0\")\n",
    "    if not yscale:\n",
    "        ax.set_ylim(*ylim)\n",
    "    ax.set_xlabel(\"$^{15}$N (ppm)\")\n",
    "    ax.set_ylabel(\"Intensity (a.u.)\")\n",
    "    fig.canvas.draw()\n",
    "\n",
    "\n",
    "# create control elements\n",
    "yscale_slider = widgets.Checkbox(value=True, description=\"Rescale\")\n",
    "b0_slider = widgets.FloatSlider(\n",
    "    value=18.7, min=0.0, max=30.0, step=0.1, description=\"B0 [T]\"\n",
    ")\n",
    "kex_slider = widgets.FloatLogSlider(\n",
    "    value=200, min=0, max=6, step=0.01, description=\"kex [/s]\"\n",
    ")\n",
    "pb_slider = widgets.FloatSlider(\n",
    "    value=0.25, min=0.0, max=1.0, step=0.0001, description=\"pb\"\n",
    ")\n",
    "δa_slider = widgets.FloatSlider(\n",
    "    value=-2.0, min=-7.5, max=7.5, step=0.1, description=\"δa [ppm]\"\n",
    ")\n",
    "δb_slider = widgets.FloatSlider(\n",
    "    value=2.0, min=-7.5, max=7.5, step=0.1, description=\"δb [ppm]\"\n",
    ")\n",
    "r2a_slider = widgets.FloatSlider(\n",
    "    value=10.0, min=0.0, max=100.0, step=0.1, description=\"R2a [/s]\"\n",
    ")\n",
    "r2b_slider = widgets.FloatSlider(\n",
    "    value=10.0, min=0.0, max=100.0, step=0.1, description=\"R2b [/s]\"\n",
    ")\n",
    "\n",
    "kd_kon_ui = widgets.HBox([pb_slider, kex_slider])\n",
    "delta_ui = widgets.HBox([δa_slider, δb_slider])\n",
    "r2_ui = widgets.HBox([r2a_slider, r2b_slider])\n",
    "b0_l0_ui = widgets.HBox([b0_slider, yscale_slider])\n",
    "\n",
    "ui = widgets.VBox([kd_kon_ui, delta_ui, r2_ui, b0_l0_ui])\n",
    "\n",
    "out = widgets.interactive_output(\n",
    "    update_values,\n",
    "    {\n",
    "        \"yscale\": yscale_slider,\n",
    "        \"b0\": b0_slider,\n",
    "        \"kex\": kex_slider,\n",
    "        \"pb\": pb_slider,\n",
    "        \"δa\": δa_slider,\n",
    "        \"δb\": δb_slider,\n",
    "        \"r2a\": r2a_slider,\n",
    "        \"r2b\": r2b_slider,\n",
    "    },\n",
    ")\n",
    "\n",
    "display(out, ui)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "#### Case 1: Exchange Between Equally Populated States\n",
    "\n",
    "- Set the following parameters:\n",
    "  - $\\text{B}_0 = 10$ Tesla\n",
    "  - $\\delta_a = -2$ ppm\n",
    "  - $\\delta_b = 2$ ppm\n",
    "  - $R_{2a} = R_{2b} = 10$ s<sup>-1</sup>\n",
    "  - $p_b = 0.5$\n",
    "\n",
    "1. At which $\\kex$ value do we observe coalescence (the transition between slow and fast exchange)?\n",
    "2. Compare this value to the frequency difference between states **A** and **B** (the gyromagnetic ratio of nitrogen is -27.1 rad·s<sup>-1</sup>·T<sup>-1</sup>).\n",
    "3. Move $\\text{B}_0$ to 20 Tesla without changing $\\kex$. Describe the new exchange regime (slow or fast).\n",
    "4. At which $\\kex$ value do we observe coalescence (the transition between slow and fast exchange)?\n",
    "\n",
    "#### Case 2: Exchange Between Major and Minor States\n",
    "\n",
    "- Set the following parameters:\n",
    "  - $\\text{B}_0 = 10$ Tesla\n",
    "  - $\\delta_a = -2$ ppm\n",
    "  - $\\delta_b = 2$ ppm\n",
    "  - $R_{2a} = R_{2b} = 10$ s<sup>-1</sup>\n",
    "  - $p_b = 0.15$\n",
    "\n",
    "**Set $\\kex$ to 300 s<sup>-1</sup>**\n",
    "\n",
    "1. Compare the line widths of the two resonances. What are the possible consequences for the interpretation of the exchange regime?\n",
    "2. Observe how the line widths change when $\\text{B}_0$ is moved from 10 T to 20 T.\n",
    "\n",
    "**Set $\\kex$ to 3000 s<sup>-1</sup>**\n",
    "\n",
    "1. Observe how the line widths change when $\\text{B}_0$ is moved from 10 T to 20 T.\n",
    "2. Compare to the previous situation and suggest a way to distinguish between slow and fast exchange regimes in the case of highly skewed populations.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Kinetics of a Binding Equilibrium E + L $\\rightleftharpoons$ EL\n",
    "\n",
    "Describing the effect of ligand binding on the observed line shapes of a target molecule constitutes a specific case of a two-state equilibrium. The target molecule can exist in either the free or bound state, and the relative populations of these states are defined by the dissociation equilibrium constant, $\\Kd$, and the concentrations of the ligand, $L_0$, and the target molecule, $E_0$, according to the law of mass action:\n",
    "\n",
    "$$ \\Kd = \\frac{\\koff}{\\kon} = \\frac{[E][L]}{[EL]} \\tag{21} $$\n",
    "\n",
    "Here, $[E]$ and $[L]$ represent the concentrations of the free target molecule and the ligand, respectively.\n",
    "\n",
    "If we denote $\\pb$ as the bound fraction of E ($\\pb = [EL]/E_0$), $\\pb$ is one of the roots of the equation:\n",
    "\n",
    "$$\n",
    "\\pb^2 - \\pb\\left(1 + \\frac{L_0}{E_0} + \\frac{\\Kd}{E_0}\\right) + \\frac{L_0}{E_0} = 0\n",
    "\\tag{22}\n",
    "$$\n",
    "\n",
    "The $\\kex$ will be defined by the number of E molecules that form a complex per unit time and the number of complexes that dissociate at the same time:\n",
    "\n",
    "$$ \\kex = \\kon[L] + \\koff = \\kon\\left([L] + \\Kd\\right) \\tag{23} $$\n",
    "\n",
    "Equation (23) shows that the line shapes depend on the total concentration of the ligand as well as the dissociation equilibrium constants. Therefore, three cases are usually observed:\n",
    "\n",
    "### 8.1. Slow Exchange\n",
    "\n",
    "This situation occurs when the dissociation constant $\\Kd$ is very small compared to the concentration of the target molecule $E_0$ (strong binding). The formation of the complex is limited by a very slow dissociation rate, namely $\\kon \\cdot \\Kd$. During acquisition, molecules are \"frozen\" in either their bound or free states, and the two corresponding lines are visible.\n",
    "\n",
    "### 8.2. Fast Exchange\n",
    "\n",
    "When the dissociation constant $\\Kd$ is larger than the concentration of the target molecule (weak binding), the dissociation rate becomes very fast compared to the frequency difference between the free state ($\\nuf$) and the bound state ($\\nub$), leading to a single resonance whose frequency is the weighted average of these frequencies. Thus, the measurement of the frequency ($\\num$) as a function of the total ligand concentration $L_0$ provides an accurate and precise measurement of the bound fraction $\\pb$:\n",
    "\n",
    "$$ \\pb(L_0) = \\frac{\\num - \\nuf}{\\nub - \\nuf} \\tag{24} $$\n",
    "\n",
    "Fitting the experimental values of $\\pb$ as a function of $L_0$ to equation (22) provides a way to determine the dissociation equilibrium constant. Note that while $\\nuf$ is easy to determine, $\\nub$ may be challenging to measure accurately since 100% of the complex may not be reachable due to weak binding. Therefore, the non-linear fitting procedure often uses both $\\Kd$ and $\\nub$ as adjustable parameters.\n",
    "\n",
    "### 8.3. Intermediate Exchange\n",
    "\n",
    "In this situation, the on and off exchange rates match the frequency difference between the bound and free states of the observed molecule. This results in a situation where resonances vanish upon the addition of the ligand. In such cases, one can attempt to escape this regime by:\n",
    "\n",
    "- Changing the temperature (if possible).\n",
    "- Changing the static magnetic field strength, $\\text{B}_0$, or the observed nuclei to obtain a different $\\Delta \\omega$ value and subsequently define another time window.\n",
    "\n",
    "## 9. Perform Simulation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "\n",
    "\n",
    "def calculate_el(kd, e_tot, l_tot):\n",
    "    \"\"\"Computes the bound fraction\"\"\"\n",
    "    sum_ = e_tot + l_tot + kd\n",
    "    return 0.5 * (sum_ - np.sqrt(sum_**2 - 4 * e_tot * l_tot))\n",
    "\n",
    "\n",
    "# Concentration in μM\n",
    "unit = 1e-6\n",
    "\n",
    "# Fixed parameters of the titration\n",
    "e_tot = 100.0 * unit  # Concentration of protein (observed by NMR)\n",
    "\n",
    "# Compute an array of Ligand/Target ratios\n",
    "ratios = np.array([0.0, 0.25, 0.5, 0.75, 1.0, 1.5, 2]) / 2.0\n",
    "\n",
    "# NMR parameters\n",
    "spectral_width2 = 8.0  # [ppm]\n",
    "\n",
    "# set up plot\n",
    "cmap = mpl.cm.cool\n",
    "\n",
    "fig, axd = plt.subplot_mosaic(\n",
    "    [\n",
    "        [\"spectra\", \"spectra\", \"fraction\"],\n",
    "        [\"spectra\", \"spectra\", \"legend\"],\n",
    "        [\"spectra\", \"spectra\", \"legend\"],\n",
    "    ],\n",
    "    layout=\"constrained\",\n",
    "    figsize=(12, 6),\n",
    ")\n",
    "\n",
    "axd[\"legend\"].axis(\"off\")\n",
    "\n",
    "\n",
    "def update_values(b0, lmax_um, kd_um, kon, δfree, δbound, r2free, r2bound):\n",
    "    kd = kd_um * unit\n",
    "\n",
    "    # Remove old line from plot and plot new one\n",
    "    axd[\"spectra\"].cla()\n",
    "    axd[\"fraction\"].cla()\n",
    "\n",
    "    pb_list = []\n",
    "\n",
    "    for index, ratio in enumerate(ratios):\n",
    "        l_tot_um = ratio * lmax_um\n",
    "        l_tot = l_tot_um * unit\n",
    "\n",
    "        e_bound = calculate_el(kd, e_tot, l_tot)\n",
    "\n",
    "        koff = kd * kon\n",
    "        l_free = l_tot - e_bound\n",
    "        kex = kon * l_free + koff\n",
    "        pb = e_bound / e_tot\n",
    "\n",
    "        params = Parameters(b0, pb, kex, δfree, δbound, r2free, r2bound)\n",
    "        coordinates = calculate_spectrum(spectral_width2, params)\n",
    "\n",
    "        axd[\"spectra\"].plot(\n",
    "            coordinates.x,\n",
    "            coordinates.y + 0.01 * index,\n",
    "            label=f\"$L_{0}$ = {l_tot_um:6.1f} μM\",\n",
    "            color=cmap(ratio),\n",
    "            zorder=100000 - index * 10,\n",
    "        )\n",
    "        pb_list.append(pb)\n",
    "\n",
    "    axd[\"spectra\"].set_xlabel(\"$^{15}$N [ppm]\")\n",
    "    axd[\"spectra\"].set_ylabel(\"Intensity [a.u.]\")\n",
    "\n",
    "    axd[\"fraction\"].scatter(ratios * lmax_um, pb_list, c=ratios, cmap=cmap)\n",
    "    axd[\"fraction\"].set_xlabel(\"$L_{0}$ [μM]\")\n",
    "    axd[\"fraction\"].set_ylabel(\"Bound fraction\")\n",
    "    axd[\"fraction\"].set_ylim(-0.1, 1.1)\n",
    "    axd[\"fraction\"].grid(True)\n",
    "\n",
    "    axd[\"legend\"].legend(*axd[\"spectra\"].get_legend_handles_labels(), loc=\"lower left\")\n",
    "\n",
    "    fig.canvas.draw()\n",
    "\n",
    "\n",
    "# create control elements\n",
    "b0_slider = widgets.FloatSlider(value=18.7, min=0.0, max=30.0, step=0.1)\n",
    "lmax_slider = widgets.FloatLogSlider(value=2e2, min=0, max=3, step=0.01)\n",
    "kd_slider = widgets.FloatLogSlider(value=1e1, min=-3, max=3, step=0.01)\n",
    "kon_slider = widgets.FloatLogSlider(value=1e8, min=6, max=8, step=0.01)\n",
    "δfree_slider = widgets.FloatSlider(value=-2.0, min=-7.5, max=7.5, step=0.1)\n",
    "δbound_slider = widgets.FloatSlider(value=2.0, min=-7.5, max=7.5, step=0.1)\n",
    "r2free_slider = widgets.FloatSlider(value=10.0, min=0.0, max=100.0, step=0.1)\n",
    "r2bound_slider = widgets.FloatSlider(value=10.0, min=0.0, max=100.0, step=0.1)\n",
    "\n",
    "\n",
    "grid = widgets.GridspecLayout(4, 4)\n",
    "grid[0, 0] = widgets.Label(\"B0 [T]\")\n",
    "grid[0, 1] = b0_slider\n",
    "grid[0, 2] = widgets.Label(\"Lmax [μM]\")\n",
    "grid[0, 3] = lmax_slider\n",
    "grid[1, 0] = widgets.Label(\"Kd [μM]\")\n",
    "grid[1, 1] = kd_slider\n",
    "grid[1, 2] = widgets.Label(\"kon [/s]\")\n",
    "grid[1, 3] = kon_slider\n",
    "grid[2, 0] = widgets.Label(\"δ(free) [ppm]\")\n",
    "grid[2, 1] = δfree_slider\n",
    "grid[2, 2] = widgets.Label(\"δ(bound) [ppm]\")\n",
    "grid[2, 3] = δbound_slider\n",
    "grid[3, 0] = widgets.Label(\"R2(free) [/s]\")\n",
    "grid[3, 1] = r2free_slider\n",
    "grid[3, 2] = widgets.Label(\"R2(bound) [/s]\")\n",
    "grid[3, 3] = r2bound_slider\n",
    "\n",
    "out = widgets.interactive_output(\n",
    "    update_values,\n",
    "    {\n",
    "        \"b0\": b0_slider,\n",
    "        \"lmax_um\": lmax_slider,\n",
    "        \"kd_um\": kd_slider,\n",
    "        \"kon\": kon_slider,\n",
    "        \"δfree\": δfree_slider,\n",
    "        \"δbound\": δbound_slider,\n",
    "        \"r2free\": r2free_slider,\n",
    "        \"r2bound\": r2bound_slider,\n",
    "    },\n",
    ")\n",
    "\n",
    "display(grid, out)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "- Given the following study conditions, define the values of $\\Kd$ for which the exchange equilibrium is in fast exchange for all titration points:\n",
    "\n",
    "  - $\\text{B}_0 = 10$ T\n",
    "  - $R_{2,\\text{free}} = R_{2,\\text{bound}} = 10$ s<sup>-1</sup>\n",
    "  - $E_0 = 100$ µM\n",
    "\n",
    "- Comment on the initial value of $\\Delta \\delta$ (4 ppm). Is that realistic? What is the range of chemical shift perturbations usually observed? Is there a difference between proton and nitrogen chemical shifts?\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Format de la Cellule Texte Brut",
  "interpreter": {
   "hash": "ad67b4a488657d258ffafb52a3e7629a2e7cc3932451c9d5dcf159c6a51c6373"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}