diff --git a/SSBtoolkit-Tutorial1-EBRAINS.ipynb b/SSBtoolkit-Tutorial1-EBRAINS.ipynb index f168044cd983f5ed297bf310989ff981a289b0d3..35d36ef4c06e50bd0fae714746407f4f473274b5 100644 --- a/SSBtoolkit-Tutorial1-EBRAINS.ipynb +++ b/SSBtoolkit-Tutorial1-EBRAINS.ipynb @@ -24,18 +24,36 @@ }, { "cell_type": "markdown", - "id": "f0984b4a", + "id": "bc7afa75", "metadata": {}, "source": [ - "## Install Dependencies" + "## Install SSB toolkit on EBRAINS\n", + "#### If you didn't clone the repository into ebrains run the following code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5f3395d", + "metadata": {}, + "outputs": [], + "source": [ + "#@title Install dependencies\n", + "%%bash\n", + "\n", + "# download SSB source code and install dependencies\n", + "if [ ! -d \"SSBtoolkit/\" ]; then\n", + " git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet\n", + " pip install -r SSBtoolkit/requirements.txt\n", + "fi" ] }, { "cell_type": "markdown", - "id": "bc7afa75", + "id": "47f373a0", "metadata": {}, "source": [ - "### Install on SSB toolkit on Google Colab" + "#### If you have cloned the repository run the following code" ] }, { @@ -45,7 +63,7 @@ "metadata": {}, "outputs": [], "source": [ - "pip install -r SSBtoolkit/requirements.txt" + "pip install -r requirements.txt" ] }, { @@ -393,35 +411,6 @@ "💡 The potency predicted values can be exported as a python dictionary using the function `sim.PotencyToDict()` or saved in a csv file: `sim.PotencyToCSV()`. \n" ] }, - { - "cell_type": "markdown", - "id": "e9076090", - "metadata": {}, - "source": [ - "## Save results on Google Drive" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "00a87ac7", - "metadata": {}, - "outputs": [], - "source": [ - "from google.colab import drive\n", - "drive.mount('/gdrive')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6cea1600", - "metadata": {}, - "outputs": [], - "source": [ - "sim.PotencyToCSV(path='/gdrive/MyDrive/XX') ## change XX accordingly. " - ] - }, { "cell_type": "markdown", "id": "ee42b0a5", diff --git a/SSBtoolkit-Tutorial2-EBRAINS.ipynb b/SSBtoolkit-Tutorial2-EBRAINS.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..878edbdc6b6dd3f49e945029f30c50f0d2308044 --- /dev/null +++ b/SSBtoolkit-Tutorial2-EBRAINS.ipynb @@ -0,0 +1,526 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "47f44294", + "metadata": {}, + "source": [ + "<div style=\"padding-bottom:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637335206/logos/Logo_des_Forschungszentrums_J_C3_BClich_seit_2018_hcliq4.svg\" width=250 align='left' style=\"margin-top:30px\"/>\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637333672/logos/ebrains_logo_tndl7p.svg\" width=\"280\" align='left' style=\"margin-left:50px; margin-top:25px\">\n", + "</div> \n", + "\n", + "<br><br><br><br>\n", + "<br>\n", + "\n", + "# Simulation of dose-response curves of antagonists using affinity values\n", + "\n", + "In this tutorial we will simulate a mathematical model of a signaling pathway to obtain dose-response curves of antagonists, from wich their *potencies (IC$_{50}$)* can be infered. \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "## <span style=\"color:red\">This tutorial just works specifically on EBRAINS</span>\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "f3721aa6", + "metadata": {}, + "source": [ + "## Install SSB toolkit on EBRAINS\n", + "#### If you didn't clone the repository into ebrains run the following code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c240cd19", + "metadata": {}, + "outputs": [], + "source": [ + "#@title Install dependencies\n", + "%%bash\n", + "\n", + "# download SSB source code and install dependencies\n", + "if [ ! -d \"SSBtoolkit/\" ]; then\n", + " git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet\n", + " pip install -r SSBtoolkit/requirements.txt\n", + "fi" + ] + }, + { + "cell_type": "markdown", + "id": "2b70a6ff", + "metadata": {}, + "source": [ + "#### If you have cloned the repository run the following code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21685e79", + "metadata": {}, + "outputs": [], + "source": [ + "pip install -r SSBtoolkit/requirements.txt" + ] + }, + { + "cell_type": "markdown", + "id": "1d7a5524", + "metadata": {}, + "source": [ + "### Import Dependencies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b185150e", + "metadata": {}, + "outputs": [], + "source": [ + "#Import Python dependencies\n", + "import os, sys, json, math, site\n", + "import numpy as np\n", + "import pandas as pd\n", + "import seaborn as sns\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from sklearn.linear_model import LinearRegression\n", + "\n", + "from src.lib.ssbtoolkit import convert, get, binding, simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce80047b", + "metadata": {}, + "outputs": [], + "source": [ + "#Setting BioNetGen environment path\n", + "distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/'\n", + "BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:')\n", + "mypath=%env PATH\n", + "newpath=BioNetGen+mypath\n", + "%env PATH=$newpath" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "881bd7c1", + "metadata": {}, + "outputs": [], + "source": [ + "#Test bioservices\n", + "from bioservices import UniProt\n", + "u = UniProt(verbose=False)" + ] + }, + { + "cell_type": "markdown", + "id": "de2b73ec", + "metadata": {}, + "source": [ + "## Load experimental data" + ] + }, + { + "cell_type": "markdown", + "id": "f1363405", + "metadata": {}, + "source": [ + "Once the SSB environment is set up we are ready to start to simulate.\n", + "\n", + "We will begin by loading the affinity data of the Adenosine 2A receptor natural ligand (the adenosine) as well as of some antagonists. The experimental values were taken from *[Guide to Pharmacology](https://www.guidetopharmacology.org/GRAC/ObjectDisplayForward?objectId=19&familyId=3&familyType=GPCR) website*. This data can be found in `examples/A2AR_antagonists_data.csv`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04434f16", + "metadata": {}, + "outputs": [], + "source": [ + "data = pd.read_csv('example/A2AR_antagonists_data.csv')\n", + "data" + ] + }, + { + "cell_type": "markdown", + "id": "7e190d72", + "metadata": {}, + "source": [ + "In <b>Tutorial1</b> the affinity values were directly used to calculate the fraction of occupied receptors according to the occupancy theory. However, such protocol is just valid if we are dealing with agonists. In this tutorial, instead, we intend to simulate a competitive antagonism. By definition, a competitive antagonist is a drug that inhibits the action of an agonist, having no effect in the absence of the agonist.\n", + "\n", + "Therefore, applying the following equation, we can calculate the fraction of occupied receptors by the agonists in the presence of an antaogist:\n", + "\n", + "$\\frac{[L_1R]}{{R_{total}}}=\\frac{[L_1]}{[L_1]+K_{d_{L1}}\\big(1+\\frac{[L_2]}{K_{d_{L2}}}\\big)}$\n", + "\n", + "Experimentally, studies of antagonists are normally conducted through radioligand binding assays. In this kind of assays, the receptor is saturated first with the agonist, and then, aliquotes of antagonist are applied to the system in order to reproduce an inhibition curve. \n", + "\n", + "To reproduce this experiment *in silico* we have to begin with the calculation of the concentration value of the radioligand/agonist (in our case the adenosine) that saturates the receptor. To achieve this, we will use the SSBtoolkit built-in function `binding.bind()` to obtain a ligand-receptor binding curve. After we will use the `maxbend()`method to interpolate the submaximal (saturate) concentration value.\n", + "\n", + "📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f6ee67bf", + "metadata": {}, + "outputs": [], + "source": [ + "# Setting the ligand concentration range\n", + "lig_conc_min = 1E-4 # μM\n", + "lig_conc_max = 1E2 # μM\n", + "lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20)\n", + "\n", + "# Setting receptor concentration\n", + "receptor_conc = 1E-3 #μM\n", + "\n", + "# Getting adenosine pKi from experimental data\n", + "adenosine_pKi = data[data.Compound == 'adenosine'].pKi.item()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "62d2a3e1", + "metadata": {}, + "outputs": [], + "source": [ + "# To obtain the ligand-receptor binding curve for adenosine we will start by initiating a binding instance\n", + "adenosine_binding = binding()\n", + "\n", + "# Performing the binding calculation\n", + "adenosine_binding.bind(receptor_conc=receptor_conc, lig_conc_range=lig_conc_range, pKd=adenosine_pKi)\n", + "\n", + "#Extraction of the submaximal concentration value\n", + "adenosine_binding.maxbend()\n", + "\n", + "#Plot of the ligand-receptor binding curve\n", + "adenosine_binding.show_curve()" + ] + }, + { + "cell_type": "markdown", + "id": "1574b00d", + "metadata": {}, + "source": [ + "The submaximal concentration value can be accessed by the `submax_concentration` attribute" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7da98638", + "metadata": {}, + "outputs": [], + "source": [ + "submax = adenosine_binding.submax_concentration\n", + "submax" + ] + }, + { + "cell_type": "markdown", + "id": "6c08db6c", + "metadata": {}, + "source": [ + "Once we have the submaximal concentration value of adenosine that saturates the receptor we can calculate the fraction of receptors occupied by the adenosine in the presence of a range of concentration of antagonists. This calculation is automatically done during the simulation." + ] + }, + { + "cell_type": "markdown", + "id": "4cea3863", + "metadata": {}, + "source": [ + "## Selection of the signaling pathway \n", + "\n", + "The core of the SSB toolkit is, naturally, the mathematical models of the signaling pathways. \n", + "\n", + "Since G-protein sub-families are classified by their $\\alpha$ subunits, this classfication as been served to identify their signaling pathways:\n", + "* G<sub>s</sub>\n", + "* G<sub>i/o</sub>\n", + "* G<sub>q/11</sub>\n", + "* G<sub>12/13</sub>\n", + "\n", + "📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details.\n", + "\n", + "We can define manualy the G$\\alpha$ pathway we want to work with, or simply query our internal database of human GPCR pathways using the UNIPROT id of our receptor. The UNIPROT id for the human Adenosine A2 receptot is [P29274](https://www.uniprot.org/uniprot/P29274)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4082b7bc", + "metadata": {}, + "outputs": [], + "source": [ + "uniprotID = 'P29274'\n", + "gprotein = get.gprotein(uniprotID)\n", + "gprotein" + ] + }, + { + "cell_type": "markdown", + "id": "22e1a045", + "metadata": {}, + "source": [ + "<span style=\"color:black\">âš ï¸ WARNING: our toolkit was specifically design for human GPCRs. The quering for pathways for GPCRs other than Human may be included in the future. However you want to study a non-human GPCR you can still use the SSB toolkit by setting manually the G$\\alpha$ pathway.</span>" + ] + }, + { + "cell_type": "markdown", + "id": "132e0f01", + "metadata": {}, + "source": [ + "## Preparation, Simulation and Analysis\n", + "\n", + "To obtain a dose-response curve from the simulation of signaling pathways, individual simulations of the pathway according to a series of ligand concentrations must be performed (as it would be done in the wet-lab). \n", + "\n", + "To define an array of ligand concentrations we will use a geometric progression. The reason why we use a geometric progression is due to the impact of the dilution fraction on the accuracy of K<sub>d</sub> and EC<sub>50</sub>/IC<sub>50</sub> values experimentally estimated. This factor, that defines the spacing between the adjacent concentration values, has an impact on the concentration values that are on the linear portion of the curve. Therefore, using a geometric progression we can mimic the experimental conditions where each concentration equals to the power of 2 of the previous lowest concentration *([Sebaugh, J.L., 2011](https://doi.org/10.1002/pst.426))*\n", + "\n", + "<span style=\"color:black\"> âš ï¸ WARNING: the SSB toolkit uses μM as default concentration units.</span>\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc61e7c6", + "metadata": {}, + "outputs": [], + "source": [ + "# the range of ligand concentration\n", + "lig_conc_min = 1E-4 # μM\n", + "lig_conc_max = 1E4 # μM\n", + "lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values\n", + "\n", + "# Setting receptor concentration\n", + "receptor_conc = 1E-3 #μM\n", + "\n", + "# defining other simulation parameters:\n", + "time = 10000 # time of simulation in seconds\n", + "nsteps = 1000 # number of time steps" + ] + }, + { + "cell_type": "markdown", + "id": "53e52b1b", + "metadata": {}, + "source": [ + "## Integration of ODEs \n", + "\n", + "After having defined all the simulation parameters we are ready to proceed with the simulations. A simulation of a methamatical model of a signaling pathway consists of the integration of a set of ordinary differential equations (ODEs) as function of time. Since each ODE represents a molecular event of the signaling pathway, when integrated over time, such equations can describe how the concentration of species inside the model changes over time. The key point of this tutorial is the use of the drug-receptor affinity value (K<sub>d</sub>) to fire up the model. With the K<sub>d</sub> values one can calculate the fraction of receptors that are occupied by the ligand in the equilibrium and, according to the *occupancy theory*, the fraction of occupied receptors represents the concentration of activated receptors in the equilibrium *([Kenakin T., 2004 ](https://doi.org/10.1016/j.tips.2004.02.012))*. 📖 Read the [Docs](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) for more details.\n", + "\n", + "\n", + "In this tutorial we want to simulate dose-response curves of antaoginsts towards Adenosine A2 receptor, so, we will use the SSBtoolkit built-in class `simulation.inhibition()`. \n", + "\n", + "\n", + "<span style='color:black'>â„¹ï¸ If you want study agonist only follow the tutorial [Simulation of dose-responses curves from affinity values](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/SSBtoolkit-Tutorial1.ipynb).</span>\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a8c6d0e1", + "metadata": {}, + "outputs": [], + "source": [ + "# we will start by creating a simulation instance.\n", + "sim = simulation.inhibition()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "91624648", + "metadata": {}, + "outputs": [], + "source": [ + "# configuration of simulation parameters\n", + "sim.SetSimulationParameters(agonist='adenosine',\n", + " agonist_affinity=adenosine_pKi,\n", + " agonist_submaximal_conc=submax,\n", + " antagonists=data.Compound.to_list()[1:3],\n", + " antagonists_affinities=data.pKi.to_list()[1:3], \n", + " lig_conc_range=lig_conc_range, \n", + " receptor_conc=receptor_conc, \n", + " pathway=gprotein, \n", + " ttotal=time, \n", + " nsteps=nsteps, \n", + " binding_kinetics=False, \n", + " binding_kinetic_parameters=None)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b52e53c4", + "metadata": {}, + "outputs": [], + "source": [ + "#Start the simulation\n", + "sim.Run()" + ] + }, + { + "cell_type": "markdown", + "id": "2103c49c", + "metadata": {}, + "source": [ + "In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance.\n", + "\n", + "The response of a signaling pathway is, naturally, represented by the increase or decrease of one of the species described by the model. So, to predict the dose-response curve we need, firstly, to extract the maximum concentration value orbserved for one specie from each individual simulation (from the series of simulations for each ligand concentration). Then, such values will be fitted to a logistic regression. \n", + "\n", + "To achieve this, we will use the analysis function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e8dff1c8", + "metadata": {}, + "outputs": [], + "source": [ + "sim.Analysis()" + ] + }, + { + "cell_type": "markdown", + "id": "fbdd0773", + "metadata": {}, + "source": [ + "We can now to plot the dose-response curves:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4cbb0861", + "metadata": {}, + "outputs": [], + "source": [ + "sim.Curve()" + ] + }, + { + "cell_type": "markdown", + "id": "c3392404", + "metadata": {}, + "source": [ + "Finnaly, from the dose-response curves we can interpolate the IC<sub>50</sub> values.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f5f55d6d", + "metadata": {}, + "outputs": [], + "source": [ + "sim.Potency()" + ] + }, + { + "cell_type": "markdown", + "id": "e8da4874", + "metadata": {}, + "source": [ + "💡 The potency predicted values can be exported as a python dictionary using the function `sim.PotencyToDict()` or saved in a csv file: `sim.PotencyToCSV()`. " + ] + }, + { + "cell_type": "markdown", + "id": "7650013b", + "metadata": {}, + "source": [ + "## Save results on Google Drive" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2d5ed71a", + "metadata": {}, + "outputs": [], + "source": [ + "from google.colab import drive\n", + "drive.mount('/gdrive')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16081378", + "metadata": {}, + "outputs": [], + "source": [ + "sim.PotencyToCSV(path='/gdrive/MyDrive/XX') ## change XX accordingly. " + ] + }, + { + "cell_type": "markdown", + "id": "f2f97975", + "metadata": {}, + "source": [ + "## Congratulations!\n", + "\n", + "Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with SSB toolkit, we encourage you to finish the rest of the tutorials in this series. \n", + "\n", + "## Cite Us\n", + "\n", + "If you use or adapt this tutorial for your own research projects please cite us.\n", + "\n", + "```\n", + "@article{XXX,\n", + " title={XXX},\n", + " author={XXX},\n", + " publisher={XXX},\n", + " note={\\url{XXX}},\n", + " year={XXX}\n", + "}\n", + "```\n", + "\n", + "\n", + "\n", + "## Acknowledgments\n", + "\n", + "EU Human Brain Project (SGA1 and SGA2): This open source software was developed in part in the Human Brain Project funded from the European Union's Horizon 2020 Framework Programme for Research and Innovation under Specific Grant Agreements No 720270 and No. 78907 (Human Brain Project SGA1 and SGA2).\n", + "\n", + "<div style=\"padding-bottom:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + " <img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + "</div> " + ] + } + ], + "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.9.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/SSBtoolkit-Tutorial3A-EBRAINS.ipynb b/SSBtoolkit-Tutorial3A-EBRAINS.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..90d5f429c674c160b7c61dddd15bdc6ddc484739 --- /dev/null +++ b/SSBtoolkit-Tutorial3A-EBRAINS.ipynb @@ -0,0 +1,469 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "546e336c", + "metadata": {}, + "source": [ + "<div style=\"padding-bottom:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637335206/logos/Logo_des_Forschungszentrums_J_C3_BClich_seit_2018_hcliq4.svg\" width=250 align='left' style=\"margin-top:30px\"/>\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637333672/logos/ebrains_logo_tndl7p.svg\" width=\"280\" align='left' style=\"margin-left:50px; margin-top:25px\">\n", + "</div> \n", + "\n", + "<br><br><br><br>\n", + "<br>\n", + "\n", + "# Simulation of dose-response curves of agonists using kinetic values\n", + "\n", + "In this tutorial we will simulate mathematical model of a signaling pathway to obtain dose-response curves, and consequently, predict the *efficacy (EC$_{50}$)* of drugs. \n", + "\n", + "## <span style=\"color:red\">This tutorial just works specifically on EBRAINS</span>" + ] + }, + { + "cell_type": "markdown", + "id": "7f3d1a03", + "metadata": {}, + "source": [ + "## Install SSB toolkit on EBRAINS\n", + "#### If you didn't clone the repository into ebrains run the following code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a3f7f10", + "metadata": {}, + "outputs": [], + "source": [ + "#@title Install dependencies\n", + "%%bash\n", + "\n", + "# download SSB source code and install dependencies\n", + "if [ ! -d \"SSBtoolkit/\" ]; then\n", + " git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet\n", + " pip install -r SSBtoolkit/requirements.txt\n", + "fi\n" + ] + }, + { + "cell_type": "markdown", + "id": "1075b49e", + "metadata": {}, + "source": [ + "#### If you have cloned the repository run the following code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1cd931bf", + "metadata": {}, + "outputs": [], + "source": [ + "pip install -r requirements.txt" + ] + }, + { + "cell_type": "markdown", + "id": "94ae8051", + "metadata": {}, + "source": [ + "### Import Dependencies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1141a5ff", + "metadata": {}, + "outputs": [], + "source": [ + "#Import Python dependencies\n", + "import os, sys, json, math, site\n", + "import numpy as np\n", + "import pandas as pd\n", + "import seaborn as sns\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from sklearn.linear_model import LinearRegression\n", + "\n", + "from src.lib.ssbtoolkit import convert, get, simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a9e6b467", + "metadata": {}, + "outputs": [], + "source": [ + "#Setting BioNetGen environment path\n", + "distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/'\n", + "BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:')\n", + "mypath=%env PATH\n", + "newpath=BioNetGen+mypath\n", + "%env PATH=$newpath" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "022c94a2", + "metadata": {}, + "outputs": [], + "source": [ + "#Test bioservices\n", + "from bioservices import UniProt\n", + "u = UniProt(verbose=False)" + ] + }, + { + "cell_type": "markdown", + "id": "f77f902a", + "metadata": {}, + "source": [ + "## Load experimental Data\n", + "\n", + "Once the SSBtoolkit environment is set up we are ready to start to simulate.\n", + "\n", + "We will begin by loading the kinetic data of some ligands of the Adenosine 2A receptor. This data was taken from *([Guo, D. et al., 2017](https://pubs.acs.org/doi/10.1021/acs.chemrev.6b00025))*. This data can be found in `examples/A2AR_Kinetic_values_example.csv`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "50ccd993", + "metadata": {}, + "outputs": [], + "source": [ + "data = pd.read_csv('example/A2AR_Kinetic_values_example.csv')\n", + "data" + ] + }, + { + "cell_type": "markdown", + "id": "ca387573", + "metadata": {}, + "source": [ + "Commonly the kinetic parameters, k<sub>on</sub> and k<sub>off</sub>, are measured experimentally at low temperatures in order to slow down the reaction, so it can be measured. Therefore, we will use the SSBtoolkit built-in function `convert.KineticTempScale()` to scale the kinetic parameters to room temperature (25 °C). This is achieved using the free energy equation:\n", + "\n", + "$ \\Delta G = -RTln(K_d)$ \n", + "\n", + "$K_d = \\frac{k_{off}}{k_{on}}$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a8f7d20b", + "metadata": {}, + "outputs": [], + "source": [ + "#We will create a new dataframe with the scaled values \n", + "data_scaled = data.copy()\n", + "data_scaled['kon (1/(uM*s))'] = data.apply(lambda row: convert.KineticTempScale(row['kon (1/(uM*s))'],row['koff (1/s)'], row['T (°C)'], 25, Tu='C')[0], axis=1)\n", + "data_scaled['koff (1/s)'] = data.apply(lambda row: convert.KineticTempScale(row['kon (1/(uM*s))'],row['koff (1/s)'], row['T (°C)'], 25, Tu='C')[1], axis=1)\n", + "data_scaled['T (°C)'] = 25\n", + "data_scaled" + ] + }, + { + "cell_type": "markdown", + "id": "a6c4a005", + "metadata": {}, + "source": [ + "## Selection of the signaling pathway \n", + "\n", + "The core of the SSB framework is, naturally, the mathematical models of the GPCRs' signaling pathways. \n", + "\n", + "Since G-protein sub-families are classified by their $\\alpha$ subunits, this classfication as been served to identify their signaling pathways:\n", + "* G<sub>s</sub>\n", + "* G<sub>i/o</sub>\n", + "* G<sub>q/11</sub>\n", + "* G<sub>12/13</sub>\n", + "\n", + "📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details.\n", + "\n", + "We can define manualy the G$\\alpha$ pathway we want to work with, or simply query our internal database of human GPCR pathways using the UNIPROT id of our receptor using the SSBtoolkit built-in function `get.gprotein()`. The UNIPROT id for the human Adenosine A2 receptot is [P29274](https://www.uniprot.org/uniprot/P29274)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0fc98ab5", + "metadata": {}, + "outputs": [], + "source": [ + "uniprotID = 'P29274'\n", + "gprotein=get.gprotein(uniprotID)\n", + "gprotein" + ] + }, + { + "cell_type": "markdown", + "id": "4cc11bca", + "metadata": {}, + "source": [ + "<span style=\"color:black\">âš ï¸ WARNING: our toolkit was specifically design for human GPCRs. The quering for pathways for GPCRs other than Human may be included in the future. However, if you want to study a non-human GPCR you can still use the SSB toolkit by setting manually the G$\\alpha$ pathway.</span>" + ] + }, + { + "cell_type": "markdown", + "id": "bb8ea233", + "metadata": {}, + "source": [ + "## Preparation, Simulation and Analysis\n", + "\n", + "To obtain a dose-response curve from the simulation of signaling pathways, individual simulations of the pathway according to a series of ligand concentrations must be performed (as it would be done in the wet-lab). \n", + "\n", + "To define an array of ligand concentrations we will use a geometric progression. The reason why we use a geometric progression is due to the impact of the dilution fraction on the accuracy of K<sub>d</sub> and EC<sub>50</sub>/IC<sub>50</sub> values experimentally estimated. This factor, that defines the spacing between the adjacent concentration values, has an impact on the concentration values that are on the linear portion of the curve. Therefore, using a geometric progression we can mimic the experimental conditions where each concentration equals to the power of 2 of the previous lowest concentration *([Sebaugh, J.L., 2011](https://doi.org/10.1002/pst.426))*\n", + "\n", + "<span style=\"color:black\"> âš ï¸ WARNING: the SSB toolkit uses μM as default concentration units.</span>\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8177d2e3", + "metadata": {}, + "outputs": [], + "source": [ + "# setting the range of ligand concentration\n", + "lig_conc_min = 1E-4 # μM\n", + "lig_conc_max = 1E4 # μM\n", + "lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values\n", + "\n", + "# Setting receptor concentration\n", + "receptor_conc = 1E-3 #μM\n", + "\n", + "# defining other simulation parameters:\n", + "time = 10000 # time of simulation in seconds\n", + "nsteps = 1000 # number of time steps" + ] + }, + { + "cell_type": "markdown", + "id": "d7851691", + "metadata": {}, + "source": [ + "## Integration of ODEs \n", + "\n", + "After having defined all the simulation parameters we are ready to proceed with the simulations. A simulation of a methamatical model of a signaling pathway consists of the integration of a set of ordinary differential equations (ODEs) as function of time. Since each ODE represents a molecular event of the signaling pathway, when integrated over time, such equations can describe how the concentration of species inside the model changes over time.\n", + "\n", + "In the <b>Tutorial1</b>, the key point is using the drug-receptor affinity value (K<sub>d</sub>) to fire up the model. However, in this tutorial we will use the kinetic parameters to simulate the ligand-receptor binding.\n", + "\n", + "Because we are using kinetic values we have to set `kinetics=True` in the `simulation.activation.SetSimulationParameters()` instance. If we use option `kinetcs=True`, we also need to pass to the instance a dictionary of parameters. We can see a list of all parameters that can be changed for each pathway [here](https://github.com/rribeiro-sci/SSBtoolkit/tree/main/SI).\n", + "\n", + "For this tutorial we just need to include the following kinetic parameters:\n", + "* `RL_kon` for the ligand-receptor interaction forward parameter (μM<sup>-1</sup>s<sup>-1</sup>)\n", + "* `RL_koff` for the ligand-receptor interaction reverse parameter (s<sup>-1</sup>)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fbd348c1", + "metadata": {}, + "outputs": [], + "source": [ + "#preparing the parameters for each ligand:\n", + "kinetic_parameters=[]\n", + "for k, v in data_scaled.iterrows():\n", + " kinetic_parameters.append({'RL_kon':v['kon (1/(uM*s))'],'RL_koff':v['koff (1/s)']})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a2d1ab6", + "metadata": {}, + "outputs": [], + "source": [ + "#creating a simulation instance\n", + "sim = simulation.activation()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db152c02", + "metadata": {}, + "outputs": [], + "source": [ + "#Setting the simulation parameters\n", + "sim.SetSimulationParameters(ligands=data_scaled.cmpd.to_list()[0:2], \n", + " pathway=gprotein, \n", + " receptor_conc=receptor_conc, \n", + " lig_conc_range=lig_conc_range, \n", + " ttotal=time, \n", + " nsteps=nsteps, \n", + " binding_kinetics=True, \n", + " binding_kinetic_parameters=kinetic_parameters)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee0e4a36", + "metadata": {}, + "outputs": [], + "source": [ + "#Start the simulation\n", + "sim.Run()" + ] + }, + { + "cell_type": "markdown", + "id": "ef1021a2", + "metadata": {}, + "source": [ + "In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance.\n", + "\n", + "The response of a signaling pathway is, naturally, represented by the increase or decrease of one of the species described by the model. So, to predict the dose-response curve we need, firstly, to extract the maximum concentration value orbserved for one specie from each individual simulation (from the series of simulations for each ligand concentration). Then, such values will be fitted to a logistic regression. \n", + "\n", + "To achieve this, we will use the analysis function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7e425249", + "metadata": {}, + "outputs": [], + "source": [ + "sim.Analysis()" + ] + }, + { + "cell_type": "markdown", + "id": "4ad91b47", + "metadata": {}, + "source": [ + "We can now to plot the dose-response curves:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bbf053a8", + "metadata": {}, + "outputs": [], + "source": [ + "sim.Curve()" + ] + }, + { + "cell_type": "markdown", + "id": "36b11a2d", + "metadata": {}, + "source": [ + "Finnaly, from the dose-response curves we can interpolate the EC<sub>50</sub> values.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "da169f05", + "metadata": {}, + "outputs": [], + "source": [ + "sim.Potency()" + ] + }, + { + "cell_type": "markdown", + "id": "99c39d5d", + "metadata": {}, + "source": [ + "💡 The potency predicted values can be exported as a python dictionary using the function `sim.PotencyToDict()` or saved in a csv file: `sim.PotencyToCSV()`. " + ] + }, + { + "cell_type": "markdown", + "id": "fe6271ed", + "metadata": {}, + "source": [ + "## Save results on Google Drive" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c5f82273", + "metadata": {}, + "outputs": [], + "source": [ + "from google.colab import drive\n", + "drive.mount('/gdrive')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd888ea5", + "metadata": {}, + "outputs": [], + "source": [ + "sim.PotencyToCSV(path='/gdrive/MyDrive/XX') ## change XX accordingly." + ] + }, + { + "cell_type": "markdown", + "id": "a7c9bbc4", + "metadata": {}, + "source": [ + "## Congratulations! \n", + "\n", + "Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with SSBtoolkit, we encourage you to finish the rest of the tutorials in this series. \n", + "\n", + "## Cite Us\n", + "\n", + "If you use or adapt this tutorial for your own research projects please cite us.\n", + "\n", + "```\n", + "@article{XXX,\n", + " title={XXX},\n", + " author={XXX},\n", + " publisher={XXX},\n", + " note={\\url{XXX}},\n", + " year={XXX}\n", + "}\n", + "```\n", + "\n", + "\n", + "\n", + "## Acknowledgments\n", + "\n", + "EU Human Brain Project (SGA1 and SGA2): This open source software was developed in part in the Human Brain Project funded from the European Union's Horizon 2020 Framework Programme for Research and Innovation under Specific Grant Agreements No 720270 and No. 78907 (Human Brain Project SGA1 and SGA2).\n", + "\n", + "<div style=\"padding-bottom:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + " <img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + "</div> " + ] + } + ], + "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.9.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/SSBtoolkit-Tutorial3B-tauRAMD-EBRAINS.ipynb b/SSBtoolkit-Tutorial3B-tauRAMD-EBRAINS.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..a3131626a6d90f9990249dcd9bd9ccbe55443402 --- /dev/null +++ b/SSBtoolkit-Tutorial3B-tauRAMD-EBRAINS.ipynb @@ -0,0 +1,543 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "2f31bbc8", + "metadata": {}, + "source": [ + "<div style=\"padding-bottom:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637335206/logos/Logo_des_Forschungszentrums_J_C3_BClich_seit_2018_hcliq4.svg\" width=250 align='left' style=\"margin-top:30px\"/>\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637333672/logos/ebrains_logo_tndl7p.svg\" width=\"280\" align='left' style=\"margin-left:50px; margin-top:25px\">\n", + "</div> \n", + "<br><br><br><br>\n", + "<br>\n", + "\n", + "# Simulation of dose-response cruves of agonsits using data acquired with tauRAMD \n", + "\n", + "\n", + "In this notebook we will a simulate mathematical model of a signaling pathway to obtain dose-response curves using kinetic values obtained with the <b>tauRAMD</b> tool.\n", + "\n", + "<b>tauRAMD</b> is a tool developed by *Kokh et al.*:\n", + "<blockquote>\n", + "\n", + "Estimation of Drug-Target Residence Times by Ï„-Random Acceleration Molecular Dynamics Simulations\n", + "Daria B. Kokh, Marta Amaral, Joerg Bomke, Ulrich Grädler, Djordje Musil, Hans-Peter Buchstaller, Matthias K. Dreyer, Matthias Frech, Maryse Lowinski, Francois Vallee, Marc Bianciotto, Alexey Rak, and Rebecca C. Wade\n", + "Journal of Chemical Theory and Computation 2018 14 (7), 3859-3869\n", + "DOI: 10.1021/acs.jctc.8b00230\n", + " \n", + "<a href=\"https://doi.org/10.1021/acs.jctc.8b00230\" target=\"_bank\"> https://doi.org/10.1021/acs.jctc.8b002308</a>\n", + "</blockquote>\n", + "\n", + "\n", + "## <span style=\"color:red\">This tutorial just works specifically on EBRAINS</span>" + ] + }, + { + "cell_type": "markdown", + "id": "e9b06661", + "metadata": {}, + "source": [ + "## Install SSB toolkit on EBRAINS\n", + "#### If you didn't clone the repository into ebrains run the following code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "209b1b96", + "metadata": {}, + "outputs": [], + "source": [ + "#@title Install dependencies\n", + "%%bash\n", + "\n", + "# download SSB source code and install dependencies\n", + "if [ ! -d \"SSBtoolkit/\" ]; then\n", + " git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet\n", + " pip install -r SSBtoolkit/requirements.txt\n", + "fi" + ] + }, + { + "cell_type": "markdown", + "id": "7c277bb0", + "metadata": {}, + "source": [ + "#### If you have cloned the repository run the following code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "65179b49", + "metadata": {}, + "outputs": [], + "source": [ + "pip install -r requirements.txt" + ] + }, + { + "cell_type": "markdown", + "id": "60583fea", + "metadata": {}, + "source": [ + "### Import Dependencies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3a595339", + "metadata": {}, + "outputs": [], + "source": [ + "#Import Python dependencies\n", + "import os, sys, json, math, site\n", + "import numpy as np\n", + "import pandas as pd\n", + "import seaborn as sns\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from sklearn.linear_model import LinearRegression\n", + " \n", + "from src.lib.ssbtoolkit import convert, get, simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8f848491", + "metadata": {}, + "outputs": [], + "source": [ + "#Setting BioNetGen environment path\n", + "distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/'\n", + "BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:')\n", + "mypath=%env PATH\n", + "newpath=BioNetGen+mypath\n", + "%env PATH=$newpath" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02c21875", + "metadata": {}, + "outputs": [], + "source": [ + "#Test bioservices\n", + "from bioservices import UniProt\n", + "u = UniProt(verbose=False)" + ] + }, + { + "cell_type": "markdown", + "id": "97e02b59", + "metadata": {}, + "source": [ + "## Calculate residence times with tauRAMD" + ] + }, + { + "cell_type": "markdown", + "id": "18d52cea", + "metadata": {}, + "source": [ + "To calculate residence times with tauRAMD we have to have the RAMD simulations already performed.\n", + "\n", + "In this tutorial we will use the pre-performed RAMD simulations (4 replicas with 15 trajectories per replica) for 2 different agonists of the Adenosine A2A receptor:\n", + "- Adenosine\n", + "- Neca\n", + "\n", + "The times.dat files can been found in `examples/tauRAMD/A2A/`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "099d60e9", + "metadata": {}, + "outputs": [], + "source": [ + "# we will start by creating a simulation instance for different ligands.\n", + "\n", + "tauADN = get.tauRAMD()\n", + "tauNEC = get.tauRAMD()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ce3b55c", + "metadata": {}, + "outputs": [], + "source": [ + "# we can run the tRAMD for each ligand\n", + "tauADN.Run(prefix='example/tauRAMD/A2A/times_ADN')\n", + "tauNEC.Run(prefix='example/tauRAMD/A2A/times_NEC')" + ] + }, + { + "cell_type": "markdown", + "id": "0e8ae29b", + "metadata": {}, + "source": [ + "The representation of the bootstrapping output and statistics are accessible through the functions:\n", + "- `get.tauRAMD.plotRTstats()`\n", + "- `get.tauRAMD.plotRTdistribuitons()`\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "831b7312", + "metadata": {}, + "outputs": [], + "source": [ + "tauADN.plotRTstats()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30d342e2", + "metadata": {}, + "outputs": [], + "source": [ + "tauADN.plotRTdistribuitons()" + ] + }, + { + "cell_type": "markdown", + "id": "3c20c5a5", + "metadata": {}, + "source": [ + "The residence times from each replica can be accessed by the attibute `RTdataframe`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b9d94a3", + "metadata": {}, + "outputs": [], + "source": [ + "tauADN.RTdataframe" + ] + }, + { + "cell_type": "markdown", + "id": "5b80e4ef", + "metadata": {}, + "source": [ + "and the average residence time can be accessed by the attribute `RT`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c00e2b6b", + "metadata": {}, + "outputs": [], + "source": [ + "tauADN.RT" + ] + }, + { + "cell_type": "markdown", + "id": "9dd98b49", + "metadata": {}, + "source": [ + "The residence time can be converted in k<sub>off</sub> by:\n", + "\n", + "$k_{off}=\\frac{1}{RT}$\n", + "\n", + "Let's create a python dictonary with the calulated k<sub>off</sub> values fro each ligand:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3d0cb2d6", + "metadata": {}, + "outputs": [], + "source": [ + "ligands = ['ADN', 'NEC']\n", + "\n", + "kinetic_parameters=[{'RL_koff':1/tauADN.RT},{'RL_koff':1/tauNEC.RT}]" + ] + }, + { + "cell_type": "markdown", + "id": "61cbf1b9", + "metadata": {}, + "source": [ + "## Preparation, Simulation and Analysis\n", + "\n", + "To obtain a dose-response curve from the simulation of signaling pathways, individual simulations of the pathway according to a series of ligand concentrations must be performed (as it would be done in the wet-lab). \n", + "\n", + "To define an array of ligand concentrations we will use a geometric progression. The reason why we use a geometric progression is due to the impact of the dilution fraction on the accuracy of K<sub>d</sub> and EC<sub>50</sub>/IC<sub>50</sub> values experimentally estimated. This factor, that defines the spacing between the adjacent concentration values, has an impact on the concentration values that are on the linear portion of the curve. Therefore, using a geometric progression we can mimic the experimental conditions where each concentration equals to the power of 2 of the previous lowest concentration *([Sebaugh, J.L., 2011](https://doi.org/10.1002/pst.426))*\n", + "\n", + "<span style=\"color:black\"> âš ï¸ WARNING: the SSB toolkit uses μM as default concentration units.</span>\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d8030846", + "metadata": {}, + "outputs": [], + "source": [ + "# setting the range of ligand concentration\n", + "lig_conc_min = 1E-4 # μM\n", + "lig_conc_max = 1E4 # μM\n", + "lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values\n", + "\n", + "# Setting receptor concentration\n", + "receptor_conc = 1E-3 #μM\n", + "\n", + "# defining other simulation parameters:\n", + "time = 10000 # time of simulation in seconds\n", + "nsteps = 1000 # number of time steps\n" + ] + }, + { + "cell_type": "markdown", + "id": "46a7f3f1", + "metadata": {}, + "source": [ + "## Selection of the signaling pathway \n", + "\n", + "The core of the SSB framework is, naturally, the mathematical models of the GPCRs' signaling pathways. \n", + "\n", + "Since G-protein sub-families are classified by their $\\alpha$ subunits, this classfication as been served to identify their signaling pathways:\n", + "* G<sub>s</sub>\n", + "* G<sub>i/o</sub>\n", + "* G<sub>q/11</sub>\n", + "* G<sub>12/13</sub>\n", + "\n", + "📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details.\n", + "\n", + "We can define manualy the G$\\alpha$ pathway we want to work with, or simply query our internal database of human GPCR pathways using the UNIPROT id of our receptor using the SSBtoolkit built-in function `get.gprotein()`. The UNIPROT id for the human Adenosine A2 receptot is [P29274](https://www.uniprot.org/uniprot/P29274)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4db29d11", + "metadata": {}, + "outputs": [], + "source": [ + "# Selecting the signaling pathway \n", + "uniprotID = 'P29274' # Adenosine A2 receptor Uniprot ID\n", + "gprotein=get.gprotein(uniprotID)\n", + "gprotein\n" + ] + }, + { + "cell_type": "markdown", + "id": "9b79bb4e", + "metadata": {}, + "source": [ + "## Integration of ODEs \n", + "\n", + "After having defined all the simulation parameters we are ready to proceed with the simulations. A simulation of a methamatical model of a signaling pathway consists of the integration of a set of ordinary differential equations (ODEs) as function of time. Since each ODE represents a molecular event of the signaling pathway, when integrated over time, such equations can describe how the concentration of species inside the model changes over time.\n", + "\n", + "Because we are using kinetic values we have to set `kinetics=True` in the `simulation.activation.SetSimulationParameters()` instance. If we use option `kinetcs=True`, we also need to pass to the instance a dictionary of parameters, the one we had prepared before. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b726eb7", + "metadata": {}, + "outputs": [], + "source": [ + "#create a simulation instance.\n", + "sim = simulation.activation()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a91aeaea", + "metadata": {}, + "outputs": [], + "source": [ + "#Setting the simulation parameters\n", + "sim.SetSimulationParameters(ligands=ligands, \n", + " pathway=gprotein, \n", + " receptor_conc=receptor_conc, \n", + " lig_conc_range=lig_conc_range, \n", + " ttotal=time, \n", + " nsteps=nsteps, \n", + " binding_kinetics=True, \n", + " binding_kinetic_parameters=kinetic_parameters)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6a30d341", + "metadata": {}, + "outputs": [], + "source": [ + "#Running the simulation\n", + "sim.Run()" + ] + }, + { + "cell_type": "markdown", + "id": "39d2dd3c", + "metadata": {}, + "source": [ + "In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance.\n", + "\n", + "The response of a signaling pathway is, naturally, represented by the increase or decrease of one of the species described by the model. So, to predict the dose-response curve we need, firstly, to extract the maximum concentration value orbserved for one specie from each individual simulation (from the series of simulations for each ligand concentration). Then, such values will be fitted to a logistic regression. \n", + "\n", + "To achieve this, we will use the analysis attribute:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d9b59f02", + "metadata": {}, + "outputs": [], + "source": [ + "sim.Analysis()" + ] + }, + { + "cell_type": "markdown", + "id": "ebca39a4", + "metadata": {}, + "source": [ + "We can now to plot the dose-response curves:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "53bfafb5", + "metadata": {}, + "outputs": [], + "source": [ + "sim.Curve()" + ] + }, + { + "cell_type": "markdown", + "id": "f7a9ec3d", + "metadata": {}, + "source": [ + "Finnaly, from the dose-response curves we can interpolate the EC<sub>50</sub> values." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8300a29", + "metadata": {}, + "outputs": [], + "source": [ + "sim.Potency()" + ] + }, + { + "cell_type": "markdown", + "id": "37ce7d4e", + "metadata": {}, + "source": [ + "💡 The potency predicted values can be exported as a python dictionary using the function `sim.PotencyToDict()` or saved in a csv file: `sim.PotencyToCSV()`. " + ] + }, + { + "cell_type": "markdown", + "id": "e9e16ef7", + "metadata": {}, + "source": [ + "## Save results on Google Drive" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3ace9c6", + "metadata": {}, + "outputs": [], + "source": [ + "from google.colab import drive\n", + "drive.mount('/gdrive')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2df68659", + "metadata": {}, + "outputs": [], + "source": [ + "sim.PotencyToCSV(path='/gdrive/MyDrive/XX') ## change XX accordingly." + ] + }, + { + "cell_type": "markdown", + "id": "65febe37", + "metadata": {}, + "source": [ + "## Congratulations! \n", + "\n", + "Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with SSBtoolkit, we encourage you to finish the rest of the tutorials in this series. \n", + "\n", + "## Cite Us\n", + "\n", + "If you use or adapt this tutorial for your own research projects please cite us.\n", + "\n", + "```\n", + "@article{XXX,\n", + " title={XXX},\n", + " author={XXX},\n", + " publisher={XXX},\n", + " note={\\url{XXX}},\n", + " year={XXX}\n", + "}\n", + "```\n", + "\n", + "\n", + "\n", + "## Acknowledgments\n", + "\n", + "EU Human Brain Project (SGA1 and SGA2): This open source software was developed in part in the Human Brain Project funded from the European Union's Horizon 2020 Framework Programme for Research and Innovation under Specific Grant Agreements No 720270 and No. 78907 (Human Brain Project SGA1 and SGA2).\n", + "\n", + "<div style=\"padding-bottom:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + " <img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + "</div> " + ] + } + ], + "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.9.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/SSBtoolkit-Tutorial4-OXTR-EBRAINS.ipynb b/SSBtoolkit-Tutorial4-OXTR-EBRAINS.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..501680ded493ca44d9a258907877c6e84252b542 --- /dev/null +++ b/SSBtoolkit-Tutorial4-OXTR-EBRAINS.ipynb @@ -0,0 +1,370 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e532cf03", + "metadata": {}, + "source": [ + "<div style=\"padding-bottom:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637335206/logos/Logo_des_Forschungszentrums_J_C3_BClich_seit_2018_hcliq4.svg\" width=250 align='left' style=\"margin-top:30px\"/>\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637333672/logos/ebrains_logo_tndl7p.svg\" width=\"280\" align='left' style=\"margin-left:50px; margin-top:25px\">\n", + "</div> \n", + "<br><br><br><br>\n", + "<br>\n", + "\n", + "# Exploring SSB pathways associated to disease variants \n", + "\n", + "## The OXTR pathway as a study case\n", + "\n", + "In this notebook we will simulate mathematical model of a signaling pathway of the Oxytocin receptor in order to give a rational on experimental observation on a disease associated variance of the OXT receptor.\n", + "\n", + "The data used in this notebook is based on:\n", + "<blockquote>\n", + "Meyer, M., Jurek, B., Alfonso-Prieto, M. et al. Structure-function relationships of the \n", + "disease-linked A218T oxytocin receptor variant. Mol Psychiatry (2022). \n", + "<a href=\"https://doi.org/10.1038/s41380-021-01241-8\" target=\"_bank\"> https://doi.org/10.1038/s41380-021-01241-8</a>\n", + "</blockquote>\n", + "\n", + "\n", + "\n", + "## <span style=\"color:red\">This tutorial just works specifically on EBRAINS</span>\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "199d4c4e", + "metadata": {}, + "source": [ + "## Install SSB toolkit on EBRAINS\n", + "#### If you didn't clone the repository into ebrains run the following code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "626cb137", + "metadata": {}, + "outputs": [], + "source": [ + "#@title Install dependencies\n", + "%%bash\n", + "\n", + "# download SSB source code and install dependencies\n", + "if [ ! -d \"SSBtoolkit/\" ]; then\n", + " git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet\n", + " pip install -r SSBtoolkit/requirements.txt\n", + "fi" + ] + }, + { + "cell_type": "markdown", + "id": "42de6245", + "metadata": {}, + "source": [ + "#### If you have cloned the repository run the following code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21fadd5a", + "metadata": {}, + "outputs": [], + "source": [ + "pip install -r requirements.txt" + ] + }, + { + "cell_type": "markdown", + "id": "439d9e19", + "metadata": {}, + "source": [ + "### Import Dependencies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "276f8a5a", + "metadata": {}, + "outputs": [], + "source": [ + "#Import Python dependencies\n", + "import os, sys, json, math, site\n", + "import numpy as np\n", + "import pandas as pd\n", + "import seaborn as sns\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from sklearn.linear_model import LinearRegression\n", + " \n", + "from src.lib.ssbtoolkit import convert, get, simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "745c0327", + "metadata": {}, + "outputs": [], + "source": [ + "#Setting BioNetGen environment path\n", + "distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/'\n", + "BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:')\n", + "mypath=%env PATH\n", + "newpath=BioNetGen+mypath\n", + "%env PATH=$newpath" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "970d01a1", + "metadata": {}, + "outputs": [], + "source": [ + "#Test bioservices\n", + "from bioservices import UniProt\n", + "u = UniProt(verbose=False)" + ] + }, + { + "cell_type": "markdown", + "id": "55b54990", + "metadata": {}, + "source": [ + "## Introduction\n", + "\n", + "The neuropeptide oxytocin (OXT) regulates multiple social and emotional behaviours, such as bonding, reciprocal trust, aggression, fear and anxiety, both in animals and humans. It has been suggested as a biomarker and target for the treatment of Autism Spectrum Disorder (ASD) [(Mayer et al. 2022)](https://doi.org/10.1038/s41380-021-01241-8).\n", + "\n", + "The OXT exerts its function by binding to its target receptor (OXTR), a member of the Class A G-protein coupled receptors (GPCR) family.\n", + "\n", + "The mutation of an alanine residue at position 218 to a threonine on the OXTR was associated with a fenotype related to ASD.\n", + "\n", + "In vitro experiments have shown a difference on the dynamics of the intracellular Ca<sup>2+</sup> between the OXTR-WT and OXTR-A218T.\n", + "\n", + "\n", + "In order to give a better understand behind these difference we will use the SSB toolkit to fit a mathematical model of the OXTR signaling pathway to the observable experimental results.\n", + "\n", + "Since we want to simulate and study the dynamics of Ca<sup>2+</sup> release from the Endoplasmic reticulum and this internal Ca<sup>2+</sup> release is trigered by the G<sub>q</sub> pathway (which OXTR couples), we will use a combination of two prevously existing models: a model for the activation of the G<sub>q</sub> protein and a dynamic model of the IP<sub>3</sub>R action. See [(Mayer et al. 2022)](https://doi.org/10.1038/s41380-021-01241-8) for all the details.\n", + "\n", + "## Experimental Data\n", + "\n", + "<blockquote>In vitro evaluation of cellular Ca<sup>2+</sup> responses using Fura-2 Ca<sup>2+</sup> imaging (see figure <b>a</b> and <b>b</b>) revealed reduced basal cytosolic Ca<sup>2+</sup> levels in A218T cells compared to WT, both in the presence and absence of extracellular Ca<sup>2+</sup> (figure <b>c</b>). \n", + "\n", + "However, the amplitude of the OXT-induced Ca<sup>2+</sup> signal was higher in A218T cells compared with WT cells incubated in Ca<sup>2+</sup>-free Ringer (figure <b>d</b>). \n", + " \n", + "The area under the curve revealed a cell line-specific effect showing a higher increase in the cytosolic Ca<sup>2+</sup> concentration upon OXT stimulation in A218T compared to WT cells (figure <b>e</b>). \n", + " \n", + "The full width at half maximum, which reflects the kinetics of the OXT-induced Ca<sup>2+</sup> response, also differed significantly between the two cell lines irrespective of the bathing solution, indicating a prolonged OXT-induced Ca<sup>2+</sup> response in the A218T cells (figure <b>f</b>).\n", + "\n", + "\n", + "<div style=\"padding-bottom:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1640172296/SSBColab/OXT_results_vsvpns.png\" width=600/>\n", + "</div> \n", + "\n", + "</blockquote>\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "16404990", + "metadata": {}, + "source": [ + "## Modeling of the Ca<sup>2+</sup> concentration differences between WT and A218T OXTR\n", + "\n", + "Our mail goal is to use the <b>SSB toolkit</b> to rationalize the experimentally measured Ca<sup>2+</sup> concentrations with 100 nM OXT and in the absence of extracellular Ca<sup>2+</sup> (<b>figure e, blue bars</b>). \n", + "\n", + "Molecular modeling results had suggested that the A218T variant affects receptor activation (see [Mayer et al. 2022](https://doi.org/10.1038/s41380-021-01241-8) for all the details). This means that we could speculate that the variant may influence the kinetic parameters that rules the activation reaction of the receptor. \n", + "\n", + "Therefore we will simulate the OXTR signalling model by varying the forward kinetic constant of the reaction that describes Gq-protein binding to the receptor, which implicitly depends on receptor activation. \n", + "\n", + "To achienve this we will use the SSB toolkit built-in function `simulation.fitModel()`. This function allows us to iteractivelly simulate the model untill to obtain the a desired ratio of the concentration of a signalling specie between a model simulated with modified parameters and a model simulated with reference parameters." + ] + }, + { + "cell_type": "markdown", + "id": "8bc1a2ab", + "metadata": {}, + "source": [ + "At first place we need to set up the following parameters:\n", + "* **fitting parameters**\n", + "* **simulation parameters**\n", + "* **pathway parameters**\n", + "\n", + "\n", + "#### Fitting Parameters\n", + "\n", + "* expratio (flt) --> experimental signalling specie concentration ratio\n", + "* target_parameter (str) --> kinetic parameter to be modified\n", + "* maxiter (int) --> maximum number of iteration\n", + "* seed (flt) --> ramdom seed for scaling the modified parameter\n", + "* seed_incrementor (flt) --> seed incrementor (each iteration will increment the seed by this value)\n", + "* seed_decrementor (flt) --> seed decrementor (each iteration will decrement the seed by this value)\n", + "\n", + "\n", + "#### Simulation Parameters\n", + "\n", + "* ttotal (int) --> time of simulation in seconds\n", + "* nsteps (int) --> number of time steps\n", + "* pathway (str) --> name of the signalling pathway \n", + "* observable (str) --> signalling specie to be measured\n", + "\n", + "#### Pahtway Parameters\n", + "\n", + "* time_in (int) --> simulation time at wich the equation will start to be simulated \n", + "* time_out (int) --> simulation time at wich the equation will stop to be simulated" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b007dfa0", + "metadata": {}, + "outputs": [], + "source": [ + "# we will start by creating a simulation instance.\n", + "sim = simulation.fitModel()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4e1d1b92", + "metadata": {}, + "outputs": [], + "source": [ + "# configuration of simulation parameters\n", + "sim.SetSimulationParameters(ttotal=250, \n", + " nsteps=10000, \n", + " pathway='OXTR_pathway', \n", + " observable='Ca', \n", + " pathway_parameters={'time_in':50, \n", + " 'time_out':51})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4cc4a960", + "metadata": {}, + "outputs": [], + "source": [ + "#Start the simulation\n", + "sim.Run(expratio=1.13, \n", + " target_parameter='R_L_Gq_trimer_kf',\n", + " maxiter=20, \n", + " seed=0.2, \n", + " seed_incrementor=0.5, \n", + " seed_decrementor=0.1)" + ] + }, + { + "cell_type": "markdown", + "id": "bd653d12", + "metadata": {}, + "source": [ + "We can now plot the curves, amplitudes, area under the curve (AUC), and the full width at half maximum (FWHM):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d37084ed", + "metadata": {}, + "outputs": [], + "source": [ + "sim.plotCurves()" + ] + }, + { + "cell_type": "markdown", + "id": "5fd46fe3", + "metadata": {}, + "source": [ + "## Save results on Google Drive" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e476cad", + "metadata": {}, + "outputs": [], + "source": [ + "from google.colab import drive\n", + "drive.mount('/gdrive')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b98afd3b", + "metadata": {}, + "outputs": [], + "source": [ + "sim.plotCurves(save=True, filename='/gdrive/MyDrive/XX') ## change XX accordingly. " + ] + }, + { + "cell_type": "markdown", + "id": "ff3d7f41", + "metadata": {}, + "source": [ + "## Congratulations! \n", + "\n", + "Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with SSBtoolkit, we encourage you to finish the rest of the tutorials in this series. \n", + "\n", + "## Cite Us\n", + "\n", + "If you use or adapt this tutorial for your own research projects please cite us.\n", + "\n", + "```\n", + "@article{XXX,\n", + " title={XXX},\n", + " author={XXX},\n", + " publisher={XXX},\n", + " note={\\url{XXX}},\n", + " year={XXX}\n", + "}\n", + "```\n", + "\n", + "\n", + "\n", + "## Acknowledgments\n", + "\n", + "EU Human Brain Project (SGA1 and SGA2): This open source software was developed in part in the Human Brain Project funded from the European Union's Horizon 2020 Framework Programme for Research and Innovation under Specific Grant Agreements No 720270 and No. 78907 (Human Brain Project SGA1 and SGA2).\n", + "\n", + "<div style=\"padding-bottom:50px\">\n", + "<img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + " <img src=\"https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg\" width=\"300\" align='left' style=\"margin-left:50px\">\n", + "</div> " + ] + } + ], + "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.9.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/apidocs.rst b/docs/apidocs.rst index df237bd403dd15f286f3641243807d4a9b2fcdf0..62392d735bc63161f3797ea7911a70c9bbf00cd1 100644 --- a/docs/apidocs.rst +++ b/docs/apidocs.rst @@ -2,3 +2,6 @@ API Documentation ================= .. automodule:: ssbtoolkit :members: + +.. automodule:: utils + :members: \ No newline at end of file diff --git a/docs/index.rst b/docs/index.rst index bd597893bc427e56385e8948c563b657a63f8da9..16d2c764f6f97eee092841a192bf2abe04d28237 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -61,7 +61,7 @@ the Apache 2.0 license. Tutorial notebooks containing minimal working examples c Developed by ============ -.. image:: _static/FZJ_logo.svg +.. image:: https://res.cloudinary.com/djz27k5hg/image/upload/v1637333672/logos/Juelich_logo_100px_ecv2hy.png :width: 200 \ @@ -74,7 +74,7 @@ Funded by :width: 200 :alt: Alternative text -Citation |DOI for Citing pyGOMoDo| +Citation |DOI for Citing ssbtoolkit| ================================== pyGOMoDo is research software. If you make use of pyGOMoDo in scientific @@ -100,3 +100,7 @@ Indices and tables * :ref:`genindex` * :ref:`modindex` * :ref:`search` + + +.. |DOI for Citing ssbtoolkit| image:: https://img.shields.io/badge/DOI-10.1016%2Fj.bpj.2015.08.015-blue.svg + :target: https://github.com/rribeiro-sci/SSBtoolkit \ No newline at end of file