Superoxide_Spin-Lattice_Relaxivity.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import nmrglue as ng\n",
    "import scipy as sp\n",
    "import scipy.optimize as opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def org_T1(numTR,path,folder,t1):  \n",
    "    \n",
    "    '''\n",
    "    This function organizes raw spectroscopially acquired t1-vtr data, into a list comprising a nd-array\n",
    "    containing the signal intensities corresponding to each TR for every scan in that experiment. To be used\n",
    "    for T1 calculation\n",
    "    \n",
    "    Args:\n",
    "        numTR (int): number of repetitions in T1-acquisition\n",
    "        path (str): path to folder containing data under different conditions\n",
    "        folder (list): containing directory names (str) of Bruker raw data under different conditions\n",
    "        res (int) : Number of TR images that are acquired before the null point\n",
    "        t1 (float) : Estimated T1 value (used for estimating how many images are acquired before crossing null point)\n",
    "     '''\n",
    "    \n",
    "    t1Data=[]\n",
    "    thresh=0.69*t1\n",
    "    res=[i for i, val in enumerate(TR) if val <=thresh]\n",
    "\n",
    "    for i in range(len(folder)):\n",
    "        dirs=[d for d in os.listdir((os.path.join(path,folder[i])))]\n",
    "        amps=np.zeros((len(dirs),numTR))\n",
    "\n",
    "        for j in range(len(dirs)):\n",
    "            dic,data=ng.bruker.read(os.path.join(path,folder[i],dirs[j]))\n",
    "\n",
    "            #Data acquired as real and imaginary pairs\n",
    "            count=0\n",
    "            amp=np.zeros(int(data.shape[0]/2), dtype=np.complex)\n",
    "            for k in range(int(data.shape[0]/2)):\n",
    "                amp[k]=complex(data[count], data[count+1])\n",
    "                count=count+2\n",
    "            ampMag=abs(amp)\n",
    "            sep=int(ampMag.shape[0]/numTR)\n",
    "            peak=np.zeros((numTR,sep))\n",
    "            count=0\n",
    "            sepCount=sep\n",
    "\n",
    "            for k in range(numTR):\n",
    "                peak[k,:]=ampMag[count:count+sep]\n",
    "                count=count+sep\n",
    "                \n",
    "            for k in range(numTR):\n",
    "                amps[j,k]=np.max(peak[k,:])\n",
    "            \n",
    "            #null crossing point at 0.69*T1\n",
    "            res=[k for k, val in enumerate(TR) if val <=thresh]\n",
    "            lim=res[-1]+1\n",
    "            amps[j,0:lim]=amps[j,0:lim]*-1\n",
    "            amps[j,:]=amps[j,:]+abs(amps[j,:].min())\n",
    "\n",
    "        t1Data.append(amps)\n",
    "    return t1Data\n",
    "\n",
    "def t1r_vtr(TR, Mo,Minf, t1):\n",
    "    '''\n",
    "    3-parameter inverision-recovery equation\n",
    "    '''\n",
    "    return Mo+(Minf-Mo)*(1-(2*(np.exp(-TR/t1))))\n",
    "\n",
    "def T1_Spect(t1Data,TR):\n",
    "    '''\n",
    "    This function calculates the 1/T1 value of spectroscopically acquired t1-vtr data.\n",
    "    \n",
    "    Args:\n",
    "        t1Data (list): list comprising an nd-array contaning the signal intensities acquired at each TR for every\n",
    "        scan in that experiment (Output of org_T1).\n",
    "        TR (array): containing repetition times that spectroscopic t1-vtr data were acquired at.\n",
    "    '''\n",
    "    t1Result=[]\n",
    "    for i in range(len(folder)):\n",
    "        coeffs=np.zeros((t1Data[i].shape[0],3))\n",
    "\n",
    "        for j in range(t1Data[i].shape[0]):\n",
    "            b=[np.max(t1Data[i][j,:]),np.min(t1Data[i][j,:]),2.6]\n",
    "            coeffs[j,:],vmat= opt.curve_fit(t1r_vtr,TR,t1Data[i][j,:],b,method='lm')\n",
    "            \n",
    "        t1Result.append(1./coeffs[:,2])\n",
    "    return t1Result\n",
    "\n",
    "def calcROSConc(a0,t,pH):\n",
    "    \n",
    "    '''\n",
    "    This function calculates the pH and time-dependent concentration (M) of the superoxide free radical \n",
    "    in aqueous solution at a certain time, t, post-dissolution. Equation for calculating reaction rate\n",
    "    constant from Eq. (1) in: Kwon BG, Kim J, Kwon J. An Advanced Kinetic Method for HO2∙/O2-∙ Determination by Using \n",
    "    Terephthalate in the Aqueous Solution. Environmental Engineering Research 2012;17(4):205-210.\n",
    "    Args:\n",
    "        a0 (int): initial concentration (M) of superoxide in solution at t=0.\n",
    "        t (int):  time at which superoxide concentration is to be calculated (s)\n",
    "        pH (float): pH of solution containing superoxide (negative logarithm of hydrogen ion concentration)\n",
    "    '''\n",
    "\n",
    "    k3=8.3E5\n",
    "    k4=9.76E7\n",
    "    kh=1.6E-5\n",
    "    hp=10**-pH\n",
    "    ko=(k3+(k4*(kh/hp)))/(1+(kh/hp))**2\n",
    "    at=(ko*t+(1/a0))**-1\n",
    "\n",
    "    return at\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F_onewayResult(statistic=1.8780628619332693, pvalue=0.16582588320225702)\n"
     ]
    }
   ],
   "source": [
    "#Figure 3A\n",
    "'''\n",
    "This script produces the data for the plot and calculate the statistics for the data as shown in Figure 3A. \n",
    "It was used to assess the impact of pH on 1/T1. See Methods section and figure legend for more detail.\n",
    "'''\n",
    "\n",
    "#Analysis setup\n",
    "path=os.path.join(os.getcwd(),'data\\\\pH')\n",
    "numTR=8\n",
    "t1=2.6\n",
    "TR=np.array([0.1, 0.8, 1.5, 2, 2.5, 4, 8.5, 11])\n",
    "folder = [d for d in os.listdir(path) if os.path.isdir(os.path.join(path, d))]\n",
    "\n",
    "#Organizing data\n",
    "t1Data=org_T1(numTR,path,folder,t1)\n",
    "\n",
    "#Calculating 1/T1\n",
    "t1pH=T1_Spect(t1Data,TR)\n",
    "\n",
    "#One-way ANOVA\n",
    "\n",
    "print(sp.stats.f_oneway(t1pH[0],t1pH[1],t1pH[2],t1pH[3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F_onewayResult(statistic=50.520481312167426, pvalue=1.0313692823538624e-13)\n",
      "F_onewayResult(statistic=0.4229161534356605, pvalue=0.8290140330694207)\n"
     ]
    }
   ],
   "source": [
    "#Figure 3B\n",
    "'''\n",
    "This script produces the data for the plot and calculate the statistics for the data as shown in Figure 3B. \n",
    "It was used to assess the impact KO2 solution on 1/T1 and the antioxidant on KO2 solution and water 1/T1.\n",
    "See Methods section and figure legend for more detail.\n",
    "'''\n",
    "\n",
    "#Analysis setup\n",
    "\n",
    "path=os.path.join(os.getcwd(),'data\\\\relaxivity')\n",
    "numTR=8\n",
    "t1=2.6\n",
    "TR=np.array([0.1, 0.8, 1.5, 2, 2.5, 4, 8.5, 11])\n",
    "folder = [d for d in os.listdir(path) if os.path.isdir(os.path.join(path, d))]\n",
    "\n",
    "#Organizing data\n",
    "t1Data=org_T1(numTR,path,folder,t1)\n",
    "\n",
    "#Calculating 1/T1\n",
    "t1KO2=T1_Spect(t1Data,TR)\n",
    "\n",
    "#One-way ANOVA of Solutions of 1/T1 KO2 in absence of the antioxidat SOD\n",
    "print(sp.stats.f_oneway(t1KO2[0],t1KO2[1],t1KO2[2],t1KO2[3],t1KO2[4],t1KO2[5]))\n",
    "\n",
    "#One-way ANOVA of Solutions of 1/T1 KO2 in presence of the antioxidat SOD\n",
    "print(sp.stats.f_oneway(t1KO2[6],t1KO2[7],t1KO2[8],t1KO2[9],t1KO2[10],t1KO2[11]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "relaxivity is 0.2897063369312174\n",
      "interccept is -0.01438674313201152\n",
      "r^2 is 0.952663532168576\n",
      "Adjusted r^2 is 0.9428697802034538\n",
      "Standard Error is 0.03228911584076144\n",
      "p-value is 0.0008539011152287562\n"
     ]
    }
   ],
   "source": [
    "#Figure 4\n",
    "'''\n",
    "This script produces the data for the plot and calculate the statistics for the data as shown in Figure 4. \n",
    "It calculates the r1 from the data in figure 3B. The 1/T1 data of KO2 in the presence of the antioxidant SOD is subtracted\n",
    "from KO2 in the absence of the antioxidant SOD. The r1 is calculated from a straight line-curve fit of this subtraction.\n",
    "See Methods section and figure legend for more detail.\n",
    "'''\n",
    "\n",
    "#Subtracting 1/T1 data of SOD- and SOD+ KO2 data, to calculate r1 from subtraction.\n",
    "t1Sub=np.zeros((int(((len(folder))/2)),len(t1Sub)))\n",
    "for i in range(int((len(t1KO2))/2)):\n",
    "    t1Sub[i,:]=t1KO2[i]-t1KO2[i+6]\n",
    "\n",
    "t1SubMean=np.mean(t1Sub,1)\n",
    "\n",
    "'''\n",
    "Setting up and carrying out calculation of time and pH dependent\n",
    "superoxide in KO2 solutions at mid-point of 2 minute acquisition \n",
    "following 10 minutesample setup\n",
    "'''\n",
    "a0=np.linspace(4,20,5)\n",
    "#pH of 4,8,12,16 and 20 mM KO2 solutions respectively (see Table 1)\n",
    "pH=[11.61,11.91,12.04,12.15,12.43]\n",
    "t=660\n",
    "at=np.zeros(len(a0)+1)\n",
    "for i in range(len(a0)):\n",
    "    at[i+1]=calcROSConc(a0[i],t,pH[i])*1000\n",
    "    \n",
    "\n",
    "\n",
    "r1, c, r, p, stdErr = stats.linregress(at, t1SubMean)\n",
    "adjrSq=1-(((1-r**2)*(36-1))/(36-6-1))\n",
    "\n",
    "print(\"relaxivity is {}\".format(r1))\n",
    "print(\"interccept is {}\".format(c))\n",
    "print(\"r^2 is {}\".format(r**2))\n",
    "print(\"Adjusted r^2 is {}\".format(adjrSq))\n",
    "print(\"Standard Error is {}\".format(stdErr))\n",
    "print(\"p-value is {}\".format(p))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "relaxivity is 0.29634251321290717\n",
      "interccept is 0.35133374391281047\n",
      "r^2 is 0.9447393007535201\n",
      "Standard Error is 0.035835722174214095\n",
      "Adjusted r^2 is 0.9333060526335587\n",
      "p-value is 0.0011669292020169773\n"
     ]
    }
   ],
   "source": [
    "#Supporting Figure S3 (SOD-)\n",
    "\n",
    "'''\n",
    "This script produces the data for the plot and calculates the statistics for the data in Supporting Figure S3. \n",
    "Specifically, it plots the 1/T1 of the KO2 data in the absence of the antioxidant SOD and assesses whether there is a \n",
    "linear relationship between 1/T1 and KO2 cocentration in the SOD- group.\n",
    "\n",
    "'''\n",
    "\n",
    "t1KO2Mean=np.zeros((int(((len(folder))/2)),(int(((len(t1KO2))/2)))))\n",
    "for i in range(int((len(t1KO2))/2)):\n",
    "    t1KO2Mean[i,:]=t1KO2[i]\n",
    "    \n",
    "t1KO2Mean=np.mean(t1KO2Mean,1)\n",
    "\n",
    "r1, c, r, p, stdErr = stats.linregress(at, t1KO2Mean)\n",
    "adjrSq=1-(((1-r**2)*(36-1))/(36-6-1))\n",
    "\n",
    "\n",
    "print(\"relaxivity is {}\".format(r1))\n",
    "print(\"interccept is {}\".format(c))\n",
    "print(\"r^2 is {}\".format(r**2))\n",
    "print(\"Standard Error is {}\".format(stdErr))\n",
    "print(\"Adjusted r^2 is {}\".format(adjrSq))\n",
    "print(\"p-value is {}\".format(p))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "relaxivity is 0.006636176281689689\n",
      "interccept is 0.36572048704482213\n",
      "r^2 is 0.36813674225197696\n",
      "Adjusted r^2 is 0.23740641306273091\n",
      "Standard Error is 0.004347052518109875\n",
      "p-value is 0.2015680490038245\n"
     ]
    }
   ],
   "source": [
    "#Supporting Figure S3 (SOD+)\n",
    "'''\n",
    "This script produces the data for the plot and calculates the statistics for the data as shown in \n",
    "Supporting Figure S3.It assesses whether there is a linear relationship between 1/T1 and KO2 concentration \n",
    "in the SOD+ group.\n",
    "\n",
    "'''\n",
    "\n",
    "t1KO2SOD=np.zeros((int(((len(folder))/2)),len(t1KO2SOD)))\n",
    "\n",
    "\n",
    "for i in range(int((len(t1KO2))/2)):\n",
    "    t1KO2SOD[i,:]=t1KO2[i+6]\n",
    "    \n",
    "t1KO2SOD=np.mean(t1KO2SOD,1)\n",
    "\n",
    "r1, c, r, p, stdErr = stats.linregress(at, t1KO2SOD)\n",
    "adjrSq=1-(((1-r**2)*(36-1))/(36-6-1))\n",
    "\n",
    "\n",
    "print(\"relaxivity is {}\".format(r1))\n",
    "print(\"interccept is {}\".format(c))\n",
    "print(\"r^2 is {}\".format(r**2))\n",
    "print(\"Adjusted r^2 is {}\".format(adjrSq))\n",
    "print(\"Standard Error is {}\".format(stdErr))\n",
    "print(\"p-value is {}\".format(p))"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}