Skip to content
Snippets Groups Projects
Commit cf87bb38 authored by Rui Ribeiro's avatar Rui Ribeiro
Browse files

release

parent d9333ca6
No related branches found
No related tags found
No related merge requests found
%% Cell type:markdown id:ebbbd7d9 tags: %% Cell type:markdown id:ebbbd7d9 tags:
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<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"/> <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"/>
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px">
<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"> <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">
</div> </div>
<br><br><br><br> <br><br><br><br>
# Simulation of dose-response curves of agonists using affinity values # Simulation of dose-response curves of agonists using affinity values
In this tutorial we will simulate a mathematical model of a signaling pathway to obtain dose-response curves, from wich *potency (EC<sub>50</sub>)* values of agonists can be infered. In this tutorial we will simulate a mathematical model of a signaling pathway to obtain dose-response curves, from wich *potency (EC<sub>50</sub>)* values of agonists can be infered.
## <span style="color:red">This tutorial just works specifically on EBRAINS</span> ## <span style="color:red">This tutorial just works specifically on EBRAINS</span>
%% Cell type:markdown id:bc7afa75 tags: %% Cell type:markdown id:bc7afa75 tags:
## Install SSB toolkit on EBRAINS ## Install SSB toolkit on EBRAINS
#### If you didn't clone the repository into ebrains run the following code #### If you didn't clone the repository into ebrains run the following code
%% Cell type:code id:e5f3395d tags: %% Cell type:code id:e5f3395d tags:
``` python ``` python
#@title Install dependencies #@title Install dependencies
%%bash %%bash
# download SSB source code and install dependencies # download SSB source code and install dependencies
if [ ! -d "SSBtoolkit/" ]; then if [ ! -d "SSBtoolkit/" ]; then
git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet
pip install -r SSBtoolkit/requirements.txt pip install -r SSBtoolkit/requirements.txt
fi fi
``` ```
%% Cell type:markdown id:47f373a0 tags: %% Cell type:markdown id:47f373a0 tags:
#### If you have cloned the repository run the following code #### If you have cloned the repository run the following code
%% Cell type:code id:d7c956eb tags: %% Cell type:code id:d7c956eb tags:
``` python ``` python
from Ipython.display import clear_ouput from Ipython.display import clear_output
pip install -r requirements.txt !pip install -r requirements.txt
clear_ouput() clear_ouput()
``` ```
%% Cell type:markdown id:0d2bec33 tags: %% Cell type:markdown id:0d2bec33 tags:
### Import Dependencies ### Import Dependencies
%% Cell type:code id:e344d8de tags: %% Cell type:code id:e344d8de tags:
``` python ``` python
#Import Python dependencies #Import Python dependencies
import os, sys, json, math, site import os, sys, json, math, site
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import seaborn as sns import seaborn as sns
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression from sklearn.linear_model import LinearRegression
from src.lib.ssbtoolkit import convert, get, simulation from src.lib.ssbtoolkit import convert, get, simulation
``` ```
%% Cell type:code id:c123810d tags: %% Cell type:code id:c123810d tags:
``` python ``` python
#Setting BioNetGen environment path #Setting BioNetGen environment path
distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/' distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/'
BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:') BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:')
mypath=%env PATH mypath=%env PATH
newpath=BioNetGen+mypath newpath=BioNetGen+mypath
%env PATH=$newpath %env PATH=$newpath
``` ```
%% Cell type:code id:66d54bde tags: %% Cell type:code id:66d54bde tags:
``` python ``` python
#Test bioservices #Test bioservices
from bioservices import UniProt from bioservices import UniProt
u = UniProt(verbose=False) u = UniProt(verbose=False)
``` ```
%% Cell type:markdown id:ec421692 tags: %% Cell type:markdown id:ec421692 tags:
## Load experimental data ## Load experimental data
%% Cell type:markdown id:184abec2 tags: %% Cell type:markdown id:184abec2 tags:
Once the SSBtoolkit environment is set up we are ready to start to simulate. Once the SSBtoolkit environment is set up we are ready to start to simulate.
We will begin by loading the affinity data of some selective agonists of the Adenosine 2A receptor. The experimental affinity data was taken from *([Varani, K. et al., 2015](https://doi.org/10.1021/acs.jmedchem.5b00215))*. This data can be found in `examples/data.csv`. We will begin by loading the affinity data of some selective agonists of the Adenosine 2A receptor. The experimental affinity data was taken from *([Varani, K. et al., 2015](https://doi.org/10.1021/acs.jmedchem.5b00215))*. This data can be found in `examples/data.csv`.
%% Cell type:code id:be92204e tags: %% Cell type:code id:be92204e tags:
``` python ``` python
data = pd.read_csv('example/A2AR_binding_data.csv') data = pd.read_csv('example/A2AR_binding_data.csv')
data data
``` ```
%% Cell type:markdown id:b75c81c0 tags: %% Cell type:markdown id:b75c81c0 tags:
Commonly, the experimental affinity values come in different "flavors". Commonly, the experimental affinity values come in different "flavors".
The SSB framework just accepts affinity values as pK<sub>d</sub>. Since the affinity values in our data set is repported as K<sub>i</sub> and EC<sub>50</sub> and we need to convert them to pK<sub>i</sub> and pEC<sub>50</sub> repectively. The SSB framework just accepts affinity values as pK<sub>d</sub>. Since the affinity values in our data set is repported as K<sub>i</sub> and EC<sub>50</sub> and we need to convert them to pK<sub>i</sub> and pEC<sub>50</sub> repectively.
%% Cell type:code id:b26bfc66 tags: %% Cell type:code id:b26bfc66 tags:
``` python ``` python
data['pKi'] = data.Ki.apply(lambda x: -np.log10(x*1E-6)) data['pKi'] = data.Ki.apply(lambda x: -np.log10(x*1E-6))
data['pEC50']=data.EC50.apply(lambda x: -np.log10(x*1E-6)) data['pEC50']=data.EC50.apply(lambda x: -np.log10(x*1E-6))
data data
``` ```
%% Cell type:markdown id:05c44c9d tags: %% Cell type:markdown id:05c44c9d tags:
## Selection of the signaling pathway ## Selection of the signaling pathway
The core of the SSB framework is, naturally, the mathematical models of the GPCRs' signaling pathways. The core of the SSB framework is, naturally, the mathematical models of the GPCRs' signaling pathways.
Since G-protein sub-families are classified by their $\alpha$ subunits, this classfication as been served to identify their signaling pathways: Since G-protein sub-families are classified by their $\alpha$ subunits, this classfication as been served to identify their signaling pathways:
* G<sub>s</sub> * G<sub>s</sub>
* G<sub>i/o</sub> * G<sub>i/o</sub>
* G<sub>q/11</sub> * G<sub>q/11</sub>
* G<sub>12/13</sub> * G<sub>12/13</sub>
📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details. 📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details.
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). 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 id:a6d7dcb1 tags: %% Cell type:code id:a6d7dcb1 tags:
``` python ``` python
uniprotID = 'P29274' uniprotID = 'P29274'
gprotein = get.gprotein(uniprotID) gprotein = get.gprotein(uniprotID)
gprotein gprotein
``` ```
%% Cell type:markdown id:9393b0cb tags: %% Cell type:markdown id:9393b0cb tags:
<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> <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:b555bc6f tags: %% Cell type:markdown id:b555bc6f tags:
## Preparation, Simulation and Analysis ## Preparation, Simulation and Analysis
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). 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).
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))* 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))*
<span style="color:black"> ⚠️ WARNING: the SSB toolkit uses μM as default concentration units.</span> <span style="color:black"> ⚠️ WARNING: the SSB toolkit uses μM as default concentration units.</span>
%% Cell type:code id:252c3fe8 tags: %% Cell type:code id:252c3fe8 tags:
``` python ``` python
# setting the range of ligand concentrations # setting the range of ligand concentrations
lig_conc_min = 1E-4 # μM lig_conc_min = 1E-4 # μM
lig_conc_max = 1E4 # μM lig_conc_max = 1E4 # μM
lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values
``` ```
%% Cell type:code id:a4cf0857 tags: %% Cell type:code id:a4cf0857 tags:
``` python ``` python
lig_conc_range lig_conc_range
``` ```
%% Cell type:markdown id:0dbc3673 tags: %% Cell type:markdown id:0dbc3673 tags:
SSB simulations are also sensible to the concentration of the protein, wich must be, again, in <b><i>μM</i></b>. However, the concentration values reported in functional assays comes, commonly, in <b><i>μg</i></b> of protein. SSB simulations are also sensible to the concentration of the protein, wich must be, again, in <b><i>μM</i></b>. However, the concentration values reported in functional assays comes, commonly, in <b><i>μg</i></b> of protein.
The experimental data we are using was obtained from a functional assay using 50 <b><i>μg</i></b> of protein. We can use the SSBtoolkit built-in function `convert.microgr2nanomolar()` to convert <b><i>μg</i></b> of protein in <b><i>μM</i></b>. The experimental data we are using was obtained from a functional assay using 50 <b><i>μg</i></b> of protein. We can use the SSBtoolkit built-in function `convert.microgr2nanomolar()` to convert <b><i>μg</i></b> of protein in <b><i>μM</i></b>.
%% Cell type:code id:ca0b4a08 tags: %% Cell type:code id:ca0b4a08 tags:
``` python ``` python
# the following function takes as input the UniProtID and the concentration in μg. # the following function takes as input the UniProtID and the concentration in μg.
prot_conc = convert.microgr2nanomolar(uniprotID, 50) prot_conc = convert.microgr2nanomolar(uniprotID, 50)
print(prot_conc, 'μM') print(prot_conc, 'μM')
``` ```
%% Cell type:markdown id:8f795fe7 tags: %% Cell type:markdown id:8f795fe7 tags:
Finnaly, we need to define the length and number of time steps (resolution) of the simulation. Finnaly, we need to define the length and number of time steps (resolution) of the simulation.
%% Cell type:code id:4cdacd1a tags: %% Cell type:code id:4cdacd1a tags:
``` python ``` python
time = 10000 # time of simulation in seconds time = 10000 # time of simulation in seconds
nsteps = 1000 # number of time steps nsteps = 1000 # number of time steps
``` ```
%% Cell type:markdown id:ba2bc317 tags: %% Cell type:markdown id:ba2bc317 tags:
## Integration of ODEs ## Integration of ODEs
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. 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.
In this tutorial we want to simulate dose-response curves of agonists towards Adenosine A2 receptor, so, we will use the SSBtoolkit built-in class `simulation.activation()`. In this tutorial we want to simulate dose-response curves of agonists towards Adenosine A2 receptor, so, we will use the SSBtoolkit built-in class `simulation.activation()`.
<span style='color:black'>ℹ️ If you want study antagonists follow the tutorial [Simulation of dose-responses curves of antagonists](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/SSBtoolkit_Tutorial2.ipynb).</span> <span style='color:black'>ℹ️ If you want study antagonists follow the tutorial [Simulation of dose-responses curves of antagonists](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/SSBtoolkit_Tutorial2.ipynb).</span>
%% Cell type:code id:5b3bf63d tags: %% Cell type:code id:5b3bf63d tags:
``` python ``` python
# we will start by creating a simulation instance. # we will start by creating a simulation instance.
sim = simulation.activation() sim = simulation.activation()
``` ```
%% Cell type:code id:c1ae3b1f tags: %% Cell type:code id:c1ae3b1f tags:
``` python ``` python
# configuration of simulation parameters # configuration of simulation parameters
sim.SetSimulationParameters(ligands=data.id.to_list()[:2], sim.SetSimulationParameters(ligands=data.id.to_list()[:2],
affinities=data.pKi.to_list()[:2], affinities=data.pKi.to_list()[:2],
pathway=gprotein, pathway=gprotein,
receptor_conc=prot_conc, receptor_conc=prot_conc,
lig_conc_range=lig_conc_range, lig_conc_range=lig_conc_range,
ttotal=time, ttotal=time,
nsteps=nsteps, nsteps=nsteps,
binding_kinetics=False) binding_kinetics=False)
``` ```
%% Cell type:code id:5f8df280 tags: %% Cell type:code id:5f8df280 tags:
``` python ``` python
#Start the simulation #Start the simulation
sim.Run() sim.Run()
``` ```
%% Cell type:markdown id:697afeef tags: %% Cell type:markdown id:697afeef tags:
In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance. In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance.
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. 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.
To achieve this, we will use the analysis function: To achieve this, we will use the analysis function:
%% Cell type:code id:8b21a9a2 tags: %% Cell type:code id:8b21a9a2 tags:
``` python ``` python
sim.Analysis() sim.Analysis()
``` ```
%% Cell type:markdown id:37fb6a47 tags: %% Cell type:markdown id:37fb6a47 tags:
We can now to plot the dose-response curves: We can now to plot the dose-response curves:
%% Cell type:code id:6aeba621 tags: %% Cell type:code id:6aeba621 tags:
``` python ``` python
sim.Curve() sim.Curve()
``` ```
%% Cell type:markdown id:5291e364 tags: %% Cell type:markdown id:5291e364 tags:
Finnaly, from the dose-response curves we can interpolate the EC<sub>50</sub> values. Finnaly, from the dose-response curves we can interpolate the EC<sub>50</sub> values.
%% Cell type:code id:dcc07180 tags: %% Cell type:code id:dcc07180 tags:
``` python ``` python
sim.Potency() sim.Potency()
``` ```
%% Cell type:markdown id:561e9437 tags: %% Cell type:markdown id:561e9437 tags:
💡 The potency predicted values can be exported as a python dictionary using the function `sim.PotencyToDict()` or saved in a csv file: `sim.PotencyToCSV()`. 💡 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:ee42b0a5 tags: %% Cell type:markdown id:ee42b0a5 tags:
## Congratulations! ## Congratulations!
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. 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.
## Cite Us ## Cite Us
If you use or adapt this tutorial for your own research projects please cite us. If you use or adapt this tutorial for your own research projects please cite us.
``` ```
@article{XXX, @article{XXX,
title={XXX}, title={XXX},
author={XXX}, author={XXX},
publisher={XXX}, publisher={XXX},
note={\url{XXX}}, note={\url{XXX}},
year={XXX} year={XXX}
} }
``` ```
## Acknowledgments ## Acknowledgments
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). 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).
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg" width="300" align='left' style="margin-left:50px">
</div> </div>
......
%% Cell type:markdown id:47f44294 tags: %% Cell type:markdown id:47f44294 tags:
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<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"/> <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"/>
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px">
<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"> <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">
</div> </div>
<br><br><br><br> <br><br><br><br>
<br> <br>
# Simulation of dose-response curves of antagonists using affinity values # Simulation of dose-response curves of antagonists using affinity values
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. 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.
## <span style="color:red">This tutorial just works specifically on EBRAINS</span> ## <span style="color:red">This tutorial just works specifically on EBRAINS</span>
%% Cell type:markdown id:f3721aa6 tags: %% Cell type:markdown id:f3721aa6 tags:
## Install SSB toolkit on EBRAINS ## Install SSB toolkit on EBRAINS
#### If you didn't clone the repository into ebrains run the following code #### If you didn't clone the repository into ebrains run the following code
%% Cell type:code id:c240cd19 tags: %% Cell type:code id:c240cd19 tags:
``` python ``` python
#@title Install dependencies #@title Install dependencies
%%bash %%bash
# download SSB source code and install dependencies # download SSB source code and install dependencies
if [ ! -d "SSBtoolkit/" ]; then if [ ! -d "SSBtoolkit/" ]; then
git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet
pip install -r SSBtoolkit/requirements.txt pip install -r SSBtoolkit/requirements.txt
fi fi
``` ```
%% Cell type:markdown id:2b70a6ff tags: %% Cell type:markdown id:2b70a6ff tags:
#### If you have cloned the repository run the following code #### If you have cloned the repository run the following code
%% Cell type:code id:21685e79 tags: %% Cell type:code id:21685e79 tags:
``` python ``` python
from Ipython.display import clear_ouput from Ipython.display import clear_output
pip install -r requirements.txt !pip install -r requirements.txt
clear_ouput() clear_ouput()
``` ```
%% Cell type:markdown id:1d7a5524 tags: %% Cell type:markdown id:1d7a5524 tags:
### Import Dependencies ### Import Dependencies
%% Cell type:code id:b185150e tags: %% Cell type:code id:b185150e tags:
``` python ``` python
#Import Python dependencies #Import Python dependencies
import os, sys, json, math, site import os, sys, json, math, site
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import seaborn as sns import seaborn as sns
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression from sklearn.linear_model import LinearRegression
from src.lib.ssbtoolkit import convert, get, binding, simulation from src.lib.ssbtoolkit import convert, get, binding, simulation
``` ```
%% Cell type:code id:ce80047b tags: %% Cell type:code id:ce80047b tags:
``` python ``` python
#Setting BioNetGen environment path #Setting BioNetGen environment path
distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/' distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/'
BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:') BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:')
mypath=%env PATH mypath=%env PATH
newpath=BioNetGen+mypath newpath=BioNetGen+mypath
%env PATH=$newpath %env PATH=$newpath
``` ```
%% Cell type:code id:881bd7c1 tags: %% Cell type:code id:881bd7c1 tags:
``` python ``` python
#Test bioservices #Test bioservices
from bioservices import UniProt from bioservices import UniProt
u = UniProt(verbose=False) u = UniProt(verbose=False)
``` ```
%% Cell type:markdown id:de2b73ec tags: %% Cell type:markdown id:de2b73ec tags:
## Load experimental data ## Load experimental data
%% Cell type:markdown id:f1363405 tags: %% Cell type:markdown id:f1363405 tags:
Once the SSB environment is set up we are ready to start to simulate. Once the SSB environment is set up we are ready to start to simulate.
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`. 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 id:04434f16 tags: %% Cell type:code id:04434f16 tags:
``` python ``` python
data = pd.read_csv('example/A2AR_antagonists_data.csv') data = pd.read_csv('example/A2AR_antagonists_data.csv')
data data
``` ```
%% Cell type:markdown id:7e190d72 tags: %% Cell type:markdown id:7e190d72 tags:
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. 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.
Therefore, applying the following equation, we can calculate the fraction of occupied receptors by the agonists in the presence of an antaogist: Therefore, applying the following equation, we can calculate the fraction of occupied receptors by the agonists in the presence of an antaogist:
$\frac{[L_1R]}{{R_{total}}}=\frac{[L_1]}{[L_1]+K_{d_{L1}}\big(1+\frac{[L_2]}{K_{d_{L2}}}\big)}$ $\frac{[L_1R]}{{R_{total}}}=\frac{[L_1]}{[L_1]+K_{d_{L1}}\big(1+\frac{[L_2]}{K_{d_{L2}}}\big)}$
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. 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.
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. 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.
📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details. 📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details.
%% Cell type:code id:f6ee67bf tags: %% Cell type:code id:f6ee67bf tags:
``` python ``` python
# Setting the ligand concentration range # Setting the ligand concentration range
lig_conc_min = 1E-4 # μM lig_conc_min = 1E-4 # μM
lig_conc_max = 1E2 # μM lig_conc_max = 1E2 # μM
lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20)
# Setting receptor concentration # Setting receptor concentration
receptor_conc = 1E-3 #μM receptor_conc = 1E-3 #μM
# Getting adenosine pKi from experimental data # Getting adenosine pKi from experimental data
adenosine_pKi = data[data.Compound == 'adenosine'].pKi.item() adenosine_pKi = data[data.Compound == 'adenosine'].pKi.item()
``` ```
%% Cell type:code id:62d2a3e1 tags: %% Cell type:code id:62d2a3e1 tags:
``` python ``` python
# To obtain the ligand-receptor binding curve for adenosine we will start by initiating a binding instance # To obtain the ligand-receptor binding curve for adenosine we will start by initiating a binding instance
adenosine_binding = binding() adenosine_binding = binding()
# Performing the binding calculation # Performing the binding calculation
adenosine_binding.bind(receptor_conc=receptor_conc, lig_conc_range=lig_conc_range, pKd=adenosine_pKi) adenosine_binding.bind(receptor_conc=receptor_conc, lig_conc_range=lig_conc_range, pKd=adenosine_pKi)
#Extraction of the submaximal concentration value #Extraction of the submaximal concentration value
adenosine_binding.maxbend() adenosine_binding.maxbend()
#Plot of the ligand-receptor binding curve #Plot of the ligand-receptor binding curve
adenosine_binding.show_curve() adenosine_binding.show_curve()
``` ```
%% Cell type:markdown id:1574b00d tags: %% Cell type:markdown id:1574b00d tags:
The submaximal concentration value can be accessed by the `submax_concentration` attribute The submaximal concentration value can be accessed by the `submax_concentration` attribute
%% Cell type:code id:7da98638 tags: %% Cell type:code id:7da98638 tags:
``` python ``` python
submax = adenosine_binding.submax_concentration submax = adenosine_binding.submax_concentration
submax submax
``` ```
%% Cell type:markdown id:6c08db6c tags: %% Cell type:markdown id:6c08db6c tags:
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. 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 tags: %% Cell type:markdown id:4cea3863 tags:
## Selection of the signaling pathway ## Selection of the signaling pathway
The core of the SSB toolkit is, naturally, the mathematical models of the signaling pathways. The core of the SSB toolkit is, naturally, the mathematical models of the signaling pathways.
Since G-protein sub-families are classified by their $\alpha$ subunits, this classfication as been served to identify their signaling pathways: Since G-protein sub-families are classified by their $\alpha$ subunits, this classfication as been served to identify their signaling pathways:
* G<sub>s</sub> * G<sub>s</sub>
* G<sub>i/o</sub> * G<sub>i/o</sub>
* G<sub>q/11</sub> * G<sub>q/11</sub>
* G<sub>12/13</sub> * G<sub>12/13</sub>
📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details. 📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details.
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). 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 id:4082b7bc tags: %% Cell type:code id:4082b7bc tags:
``` python ``` python
uniprotID = 'P29274' uniprotID = 'P29274'
gprotein = get.gprotein(uniprotID) gprotein = get.gprotein(uniprotID)
gprotein gprotein
``` ```
%% Cell type:markdown id:22e1a045 tags: %% Cell type:markdown id:22e1a045 tags:
<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> <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 tags: %% Cell type:markdown id:132e0f01 tags:
## Preparation, Simulation and Analysis ## Preparation, Simulation and Analysis
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). 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).
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))* 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))*
<span style="color:black"> ⚠️ WARNING: the SSB toolkit uses μM as default concentration units.</span> <span style="color:black"> ⚠️ WARNING: the SSB toolkit uses μM as default concentration units.</span>
%% Cell type:code id:dc61e7c6 tags: %% Cell type:code id:dc61e7c6 tags:
``` python ``` python
# the range of ligand concentration # the range of ligand concentration
lig_conc_min = 1E-4 # μM lig_conc_min = 1E-4 # μM
lig_conc_max = 1E4 # μM lig_conc_max = 1E4 # μM
lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values
# Setting receptor concentration # Setting receptor concentration
receptor_conc = 1E-3 #μM receptor_conc = 1E-3 #μM
# defining other simulation parameters: # defining other simulation parameters:
time = 10000 # time of simulation in seconds time = 10000 # time of simulation in seconds
nsteps = 1000 # number of time steps nsteps = 1000 # number of time steps
``` ```
%% Cell type:markdown id:53e52b1b tags: %% Cell type:markdown id:53e52b1b tags:
## Integration of ODEs ## Integration of ODEs
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. 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.
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()`. 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()`.
<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> <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>
%% Cell type:code id:a8c6d0e1 tags: %% Cell type:code id:a8c6d0e1 tags:
``` python ``` python
# we will start by creating a simulation instance. # we will start by creating a simulation instance.
sim = simulation.inhibition() sim = simulation.inhibition()
``` ```
%% Cell type:code id:91624648 tags: %% Cell type:code id:91624648 tags:
``` python ``` python
# configuration of simulation parameters # configuration of simulation parameters
sim.SetSimulationParameters(agonist='adenosine', sim.SetSimulationParameters(agonist='adenosine',
agonist_affinity=adenosine_pKi, agonist_affinity=adenosine_pKi,
agonist_submaximal_conc=submax, agonist_submaximal_conc=submax,
antagonists=data.Compound.to_list()[1:3], antagonists=data.Compound.to_list()[1:3],
antagonists_affinities=data.pKi.to_list()[1:3], antagonists_affinities=data.pKi.to_list()[1:3],
lig_conc_range=lig_conc_range, lig_conc_range=lig_conc_range,
receptor_conc=receptor_conc, receptor_conc=receptor_conc,
pathway=gprotein, pathway=gprotein,
ttotal=time, ttotal=time,
nsteps=nsteps, nsteps=nsteps,
binding_kinetics=False, binding_kinetics=False,
binding_kinetic_parameters=None) binding_kinetic_parameters=None)
``` ```
%% Cell type:code id:b52e53c4 tags: %% Cell type:code id:b52e53c4 tags:
``` python ``` python
#Start the simulation #Start the simulation
sim.Run() sim.Run()
``` ```
%% Cell type:markdown id:2103c49c tags: %% Cell type:markdown id:2103c49c tags:
In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance. In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance.
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. 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.
To achieve this, we will use the analysis function: To achieve this, we will use the analysis function:
%% Cell type:code id:e8dff1c8 tags: %% Cell type:code id:e8dff1c8 tags:
``` python ``` python
sim.Analysis() sim.Analysis()
``` ```
%% Cell type:markdown id:fbdd0773 tags: %% Cell type:markdown id:fbdd0773 tags:
We can now to plot the dose-response curves: We can now to plot the dose-response curves:
%% Cell type:code id:4cbb0861 tags: %% Cell type:code id:4cbb0861 tags:
``` python ``` python
sim.Curve() sim.Curve()
``` ```
%% Cell type:markdown id:c3392404 tags: %% Cell type:markdown id:c3392404 tags:
Finnaly, from the dose-response curves we can interpolate the IC<sub>50</sub> values. Finnaly, from the dose-response curves we can interpolate the IC<sub>50</sub> values.
%% Cell type:code id:f5f55d6d tags: %% Cell type:code id:f5f55d6d tags:
``` python ``` python
sim.Potency() sim.Potency()
``` ```
%% Cell type:markdown id:e8da4874 tags: %% Cell type:markdown id:e8da4874 tags:
💡 The potency predicted values can be exported as a python dictionary using the function `sim.PotencyToDict()` or saved in a csv file: `sim.PotencyToCSV()`. 💡 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 tags: %% Cell type:markdown id:7650013b tags:
## Save results on Google Drive ## Save results on Google Drive
%% Cell type:code id:2d5ed71a tags: %% Cell type:code id:2d5ed71a tags:
``` python ``` python
from google.colab import drive from google.colab import drive
drive.mount('/gdrive') drive.mount('/gdrive')
``` ```
%% Cell type:code id:16081378 tags: %% Cell type:code id:16081378 tags:
``` python ``` python
sim.PotencyToCSV(path='/gdrive/MyDrive/XX') ## change XX accordingly. sim.PotencyToCSV(path='/gdrive/MyDrive/XX') ## change XX accordingly.
``` ```
%% Cell type:markdown id:f2f97975 tags: %% Cell type:markdown id:f2f97975 tags:
## Congratulations! ## Congratulations!
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. 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.
## Cite Us ## Cite Us
If you use or adapt this tutorial for your own research projects please cite us. If you use or adapt this tutorial for your own research projects please cite us.
``` ```
@article{XXX, @article{XXX,
title={XXX}, title={XXX},
author={XXX}, author={XXX},
publisher={XXX}, publisher={XXX},
note={\url{XXX}}, note={\url{XXX}},
year={XXX} year={XXX}
} }
``` ```
## Acknowledgments ## Acknowledgments
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). 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).
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg" width="300" align='left' style="margin-left:50px">
</div> </div>
......
%% Cell type:markdown id:546e336c tags: %% Cell type:markdown id:546e336c tags:
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<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"/> <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"/>
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px">
<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"> <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">
</div> </div>
<br><br><br><br> <br><br><br><br>
<br> <br>
# Simulation of dose-response curves of agonists using kinetic values # Simulation of dose-response curves of agonists using kinetic values
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. 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.
## <span style="color:red">This tutorial just works specifically on EBRAINS</span> ## <span style="color:red">This tutorial just works specifically on EBRAINS</span>
%% Cell type:markdown id:7f3d1a03 tags: %% Cell type:markdown id:7f3d1a03 tags:
## Install SSB toolkit on EBRAINS ## Install SSB toolkit on EBRAINS
#### If you didn't clone the repository into ebrains run the following code #### If you didn't clone the repository into ebrains run the following code
%% Cell type:code id:8a3f7f10 tags: %% Cell type:code id:8a3f7f10 tags:
``` python ``` python
#@title Install dependencies #@title Install dependencies
%%bash %%bash
# download SSB source code and install dependencies # download SSB source code and install dependencies
if [ ! -d "SSBtoolkit/" ]; then if [ ! -d "SSBtoolkit/" ]; then
git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet
pip install -r SSBtoolkit/requirements.txt pip install -r SSBtoolkit/requirements.txt
fi fi
``` ```
%% Cell type:markdown id:1075b49e tags: %% Cell type:markdown id:1075b49e tags:
#### If you have cloned the repository run the following code #### If you have cloned the repository run the following code
%% Cell type:code id:1cd931bf tags: %% Cell type:code id:1cd931bf tags:
``` python ``` python
from Ipython.display import clear_ouput from Ipython.display import clear_output
pip install -r requirements.txt !pip install -r requirements.txt
clear_ouput() clear_ouput()
``` ```
%% Cell type:markdown id:94ae8051 tags: %% Cell type:markdown id:94ae8051 tags:
### Import Dependencies ### Import Dependencies
%% Cell type:code id:1141a5ff tags: %% Cell type:code id:1141a5ff tags:
``` python ``` python
#Import Python dependencies #Import Python dependencies
import os, sys, json, math, site import os, sys, json, math, site
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import seaborn as sns import seaborn as sns
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression from sklearn.linear_model import LinearRegression
from src.lib.ssbtoolkit import convert, get, simulation from src.lib.ssbtoolkit import convert, get, simulation
``` ```
%% Cell type:code id:a9e6b467 tags: %% Cell type:code id:a9e6b467 tags:
``` python ``` python
#Setting BioNetGen environment path #Setting BioNetGen environment path
distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/' distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/'
BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:') BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:')
mypath=%env PATH mypath=%env PATH
newpath=BioNetGen+mypath newpath=BioNetGen+mypath
%env PATH=$newpath %env PATH=$newpath
``` ```
%% Cell type:code id:022c94a2 tags: %% Cell type:code id:022c94a2 tags:
``` python ``` python
#Test bioservices #Test bioservices
from bioservices import UniProt from bioservices import UniProt
u = UniProt(verbose=False) u = UniProt(verbose=False)
``` ```
%% Cell type:markdown id:f77f902a tags: %% Cell type:markdown id:f77f902a tags:
## Load experimental Data ## Load experimental Data
Once the SSBtoolkit environment is set up we are ready to start to simulate. Once the SSBtoolkit environment is set up we are ready to start to simulate.
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` 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 id:50ccd993 tags: %% Cell type:code id:50ccd993 tags:
``` python ``` python
data = pd.read_csv('example/A2AR_Kinetic_values_example.csv') data = pd.read_csv('example/A2AR_Kinetic_values_example.csv')
data data
``` ```
%% Cell type:markdown id:ca387573 tags: %% Cell type:markdown id:ca387573 tags:
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: 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:
$ \Delta G = -RTln(K_d)$ $ \Delta G = -RTln(K_d)$
$K_d = \frac{k_{off}}{k_{on}}$ $K_d = \frac{k_{off}}{k_{on}}$
%% Cell type:code id:a8f7d20b tags: %% Cell type:code id:a8f7d20b tags:
``` python ``` python
#We will create a new dataframe with the scaled values #We will create a new dataframe with the scaled values
data_scaled = data.copy() data_scaled = data.copy()
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) 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)
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) 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)
data_scaled['T (°C)'] = 25 data_scaled['T (°C)'] = 25
data_scaled data_scaled
``` ```
%% Cell type:markdown id:a6c4a005 tags: %% Cell type:markdown id:a6c4a005 tags:
## Selection of the signaling pathway ## Selection of the signaling pathway
The core of the SSB framework is, naturally, the mathematical models of the GPCRs' signaling pathways. The core of the SSB framework is, naturally, the mathematical models of the GPCRs' signaling pathways.
Since G-protein sub-families are classified by their $\alpha$ subunits, this classfication as been served to identify their signaling pathways: Since G-protein sub-families are classified by their $\alpha$ subunits, this classfication as been served to identify their signaling pathways:
* G<sub>s</sub> * G<sub>s</sub>
* G<sub>i/o</sub> * G<sub>i/o</sub>
* G<sub>q/11</sub> * G<sub>q/11</sub>
* G<sub>12/13</sub> * G<sub>12/13</sub>
📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details. 📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details.
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). 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 id:0fc98ab5 tags: %% Cell type:code id:0fc98ab5 tags:
``` python ``` python
uniprotID = 'P29274' uniprotID = 'P29274'
gprotein=get.gprotein(uniprotID) gprotein=get.gprotein(uniprotID)
gprotein gprotein
``` ```
%% Cell type:markdown id:4cc11bca tags: %% Cell type:markdown id:4cc11bca tags:
<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> <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 tags: %% Cell type:markdown id:bb8ea233 tags:
## Preparation, Simulation and Analysis ## Preparation, Simulation and Analysis
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). 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).
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))* 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))*
<span style="color:black"> ⚠️ WARNING: the SSB toolkit uses μM as default concentration units.</span> <span style="color:black"> ⚠️ WARNING: the SSB toolkit uses μM as default concentration units.</span>
%% Cell type:code id:8177d2e3 tags: %% Cell type:code id:8177d2e3 tags:
``` python ``` python
# setting the range of ligand concentration # setting the range of ligand concentration
lig_conc_min = 1E-4 # μM lig_conc_min = 1E-4 # μM
lig_conc_max = 1E4 # μM lig_conc_max = 1E4 # μM
lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values
# Setting receptor concentration # Setting receptor concentration
receptor_conc = 1E-3 #μM receptor_conc = 1E-3 #μM
# defining other simulation parameters: # defining other simulation parameters:
time = 10000 # time of simulation in seconds time = 10000 # time of simulation in seconds
nsteps = 1000 # number of time steps nsteps = 1000 # number of time steps
``` ```
%% Cell type:markdown id:d7851691 tags: %% Cell type:markdown id:d7851691 tags:
## Integration of ODEs ## Integration of ODEs
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. 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.
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. 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.
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). 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).
For this tutorial we just need to include the following kinetic parameters: For this tutorial we just need to include the following kinetic parameters:
* `RL_kon` for the ligand-receptor interaction forward parameter (μM<sup>-1</sup>s<sup>-1</sup>) * `RL_kon` for the ligand-receptor interaction forward parameter (μM<sup>-1</sup>s<sup>-1</sup>)
* `RL_koff` for the ligand-receptor interaction reverse parameter (s<sup>-1</sup>) * `RL_koff` for the ligand-receptor interaction reverse parameter (s<sup>-1</sup>)
%% Cell type:code id:fbd348c1 tags: %% Cell type:code id:fbd348c1 tags:
``` python ``` python
#preparing the parameters for each ligand: #preparing the parameters for each ligand:
kinetic_parameters=[] kinetic_parameters=[]
for k, v in data_scaled.iterrows(): for k, v in data_scaled.iterrows():
kinetic_parameters.append({'RL_kon':v['kon (1/(uM*s))'],'RL_koff':v['koff (1/s)']}) kinetic_parameters.append({'RL_kon':v['kon (1/(uM*s))'],'RL_koff':v['koff (1/s)']})
``` ```
%% Cell type:code id:2a2d1ab6 tags: %% Cell type:code id:2a2d1ab6 tags:
``` python ``` python
#creating a simulation instance #creating a simulation instance
sim = simulation.activation() sim = simulation.activation()
``` ```
%% Cell type:code id:db152c02 tags: %% Cell type:code id:db152c02 tags:
``` python ``` python
#Setting the simulation parameters #Setting the simulation parameters
sim.SetSimulationParameters(ligands=data_scaled.cmpd.to_list()[0:2], sim.SetSimulationParameters(ligands=data_scaled.cmpd.to_list()[0:2],
pathway=gprotein, pathway=gprotein,
receptor_conc=receptor_conc, receptor_conc=receptor_conc,
lig_conc_range=lig_conc_range, lig_conc_range=lig_conc_range,
ttotal=time, ttotal=time,
nsteps=nsteps, nsteps=nsteps,
binding_kinetics=True, binding_kinetics=True,
binding_kinetic_parameters=kinetic_parameters) binding_kinetic_parameters=kinetic_parameters)
``` ```
%% Cell type:code id:ee0e4a36 tags: %% Cell type:code id:ee0e4a36 tags:
``` python ``` python
#Start the simulation #Start the simulation
sim.Run() sim.Run()
``` ```
%% Cell type:markdown id:ef1021a2 tags: %% Cell type:markdown id:ef1021a2 tags:
In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance. In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance.
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. 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.
To achieve this, we will use the analysis function: To achieve this, we will use the analysis function:
%% Cell type:code id:7e425249 tags: %% Cell type:code id:7e425249 tags:
``` python ``` python
sim.Analysis() sim.Analysis()
``` ```
%% Cell type:markdown id:4ad91b47 tags: %% Cell type:markdown id:4ad91b47 tags:
We can now to plot the dose-response curves: We can now to plot the dose-response curves:
%% Cell type:code id:bbf053a8 tags: %% Cell type:code id:bbf053a8 tags:
``` python ``` python
sim.Curve() sim.Curve()
``` ```
%% Cell type:markdown id:36b11a2d tags: %% Cell type:markdown id:36b11a2d tags:
Finnaly, from the dose-response curves we can interpolate the EC<sub>50</sub> values. Finnaly, from the dose-response curves we can interpolate the EC<sub>50</sub> values.
%% Cell type:code id:da169f05 tags: %% Cell type:code id:da169f05 tags:
``` python ``` python
sim.Potency() sim.Potency()
``` ```
%% Cell type:markdown id:99c39d5d tags: %% Cell type:markdown id:99c39d5d tags:
💡 The potency predicted values can be exported as a python dictionary using the function `sim.PotencyToDict()` or saved in a csv file: `sim.PotencyToCSV()`. 💡 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 tags: %% Cell type:markdown id:fe6271ed tags:
## Save results on Google Drive ## Save results on Google Drive
%% Cell type:code id:c5f82273 tags: %% Cell type:code id:c5f82273 tags:
``` python ``` python
from google.colab import drive from google.colab import drive
drive.mount('/gdrive') drive.mount('/gdrive')
``` ```
%% Cell type:code id:bd888ea5 tags: %% Cell type:code id:bd888ea5 tags:
``` python ``` python
sim.PotencyToCSV(path='/gdrive/MyDrive/XX') ## change XX accordingly. sim.PotencyToCSV(path='/gdrive/MyDrive/XX') ## change XX accordingly.
``` ```
%% Cell type:markdown id:a7c9bbc4 tags: %% Cell type:markdown id:a7c9bbc4 tags:
## Congratulations! ## Congratulations!
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. 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.
## Cite Us ## Cite Us
If you use or adapt this tutorial for your own research projects please cite us. If you use or adapt this tutorial for your own research projects please cite us.
``` ```
@article{XXX, @article{XXX,
title={XXX}, title={XXX},
author={XXX}, author={XXX},
publisher={XXX}, publisher={XXX},
note={\url{XXX}}, note={\url{XXX}},
year={XXX} year={XXX}
} }
``` ```
## Acknowledgments ## Acknowledgments
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). 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).
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg" width="300" align='left' style="margin-left:50px">
</div> </div>
......
%% Cell type:markdown id:2f31bbc8 tags: %% Cell type:markdown id:2f31bbc8 tags:
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<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"/> <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"/>
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px">
<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"> <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">
</div> </div>
<br><br><br><br> <br><br><br><br>
<br> <br>
# Simulation of dose-response cruves of agonsits using data acquired with tauRAMD # Simulation of dose-response cruves of agonsits using data acquired with tauRAMD
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. 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.
<b>tauRAMD</b> is a tool developed by *Kokh et al.*: <b>tauRAMD</b> is a tool developed by *Kokh et al.*:
<blockquote> <blockquote>
Estimation of Drug-Target Residence Times by τ-Random Acceleration Molecular Dynamics Simulations Estimation of Drug-Target Residence Times by τ-Random Acceleration Molecular Dynamics Simulations
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 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
Journal of Chemical Theory and Computation 2018 14 (7), 3859-3869 Journal of Chemical Theory and Computation 2018 14 (7), 3859-3869
DOI: 10.1021/acs.jctc.8b00230 DOI: 10.1021/acs.jctc.8b00230
<a href="https://doi.org/10.1021/acs.jctc.8b00230" target="_bank"> https://doi.org/10.1021/acs.jctc.8b002308</a> <a href="https://doi.org/10.1021/acs.jctc.8b00230" target="_bank"> https://doi.org/10.1021/acs.jctc.8b002308</a>
</blockquote> </blockquote>
## <span style="color:red">This tutorial just works specifically on EBRAINS</span> ## <span style="color:red">This tutorial just works specifically on EBRAINS</span>
%% Cell type:markdown id:e9b06661 tags: %% Cell type:markdown id:e9b06661 tags:
## Install SSB toolkit on EBRAINS ## Install SSB toolkit on EBRAINS
#### If you didn't clone the repository into ebrains run the following code #### If you didn't clone the repository into ebrains run the following code
%% Cell type:code id:209b1b96 tags: %% Cell type:code id:209b1b96 tags:
``` python ``` python
#@title Install dependencies #@title Install dependencies
%%bash %%bash
# download SSB source code and install dependencies # download SSB source code and install dependencies
if [ ! -d "SSBtoolkit/" ]; then if [ ! -d "SSBtoolkit/" ]; then
git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet
pip install -r SSBtoolkit/requirements.txt pip install -r SSBtoolkit/requirements.txt
fi fi
``` ```
%% Cell type:markdown id:7c277bb0 tags: %% Cell type:markdown id:7c277bb0 tags:
#### If you have cloned the repository run the following code #### If you have cloned the repository run the following code
%% Cell type:code id:65179b49 tags: %% Cell type:code id:65179b49 tags:
``` python ``` python
from Ipython.display import clear_ouput from Ipython.display import clear_output
pip install -r requirements.txt !pip install -r requirements.txt
clear_ouput() clear_ouput()
``` ```
%% Cell type:markdown id:60583fea tags: %% Cell type:markdown id:60583fea tags:
### Import Dependencies ### Import Dependencies
%% Cell type:code id:3a595339 tags: %% Cell type:code id:3a595339 tags:
``` python ``` python
#Import Python dependencies #Import Python dependencies
import os, sys, json, math, site import os, sys, json, math, site
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import seaborn as sns import seaborn as sns
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression from sklearn.linear_model import LinearRegression
from src.lib.ssbtoolkit import convert, get, simulation from src.lib.ssbtoolkit import convert, get, simulation
``` ```
%% Cell type:code id:8f848491 tags: %% Cell type:code id:8f848491 tags:
``` python ``` python
#Setting BioNetGen environment path #Setting BioNetGen environment path
distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/' distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/'
BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:') BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:')
mypath=%env PATH mypath=%env PATH
newpath=BioNetGen+mypath newpath=BioNetGen+mypath
%env PATH=$newpath %env PATH=$newpath
``` ```
%% Cell type:code id:02c21875 tags: %% Cell type:code id:02c21875 tags:
``` python ``` python
#Test bioservices #Test bioservices
from bioservices import UniProt from bioservices import UniProt
u = UniProt(verbose=False) u = UniProt(verbose=False)
``` ```
%% Cell type:markdown id:97e02b59 tags: %% Cell type:markdown id:97e02b59 tags:
## Calculate residence times with tauRAMD ## Calculate residence times with tauRAMD
%% Cell type:markdown id:18d52cea tags: %% Cell type:markdown id:18d52cea tags:
To calculate residence times with tauRAMD we have to have the RAMD simulations already performed. To calculate residence times with tauRAMD we have to have the RAMD simulations already performed.
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: 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:
- Adenosine - Adenosine
- Neca - Neca
The times.dat files can been found in `examples/tauRAMD/A2A/`. The times.dat files can been found in `examples/tauRAMD/A2A/`.
%% Cell type:code id:099d60e9 tags: %% Cell type:code id:099d60e9 tags:
``` python ``` python
# we will start by creating a simulation instance for different ligands. # we will start by creating a simulation instance for different ligands.
tauADN = get.tauRAMD() tauADN = get.tauRAMD()
tauNEC = get.tauRAMD() tauNEC = get.tauRAMD()
``` ```
%% Cell type:code id:7ce3b55c tags: %% Cell type:code id:7ce3b55c tags:
``` python ``` python
# we can run the tRAMD for each ligand # we can run the tRAMD for each ligand
tauADN.Run(prefix='example/tauRAMD/A2A/times_ADN') tauADN.Run(prefix='example/tauRAMD/A2A/times_ADN')
tauNEC.Run(prefix='example/tauRAMD/A2A/times_NEC') tauNEC.Run(prefix='example/tauRAMD/A2A/times_NEC')
``` ```
%% Cell type:markdown id:0e8ae29b tags: %% Cell type:markdown id:0e8ae29b tags:
The representation of the bootstrapping output and statistics are accessible through the functions: The representation of the bootstrapping output and statistics are accessible through the functions:
- `get.tauRAMD.plotRTstats()` - `get.tauRAMD.plotRTstats()`
- `get.tauRAMD.plotRTdistribuitons()` - `get.tauRAMD.plotRTdistribuitons()`
%% Cell type:code id:831b7312 tags: %% Cell type:code id:831b7312 tags:
``` python ``` python
tauADN.plotRTstats() tauADN.plotRTstats()
``` ```
%% Cell type:code id:30d342e2 tags: %% Cell type:code id:30d342e2 tags:
``` python ``` python
tauADN.plotRTdistribuitons() tauADN.plotRTdistribuitons()
``` ```
%% Cell type:markdown id:3c20c5a5 tags: %% Cell type:markdown id:3c20c5a5 tags:
The residence times from each replica can be accessed by the attibute `RTdataframe`: The residence times from each replica can be accessed by the attibute `RTdataframe`:
%% Cell type:code id:8b9d94a3 tags: %% Cell type:code id:8b9d94a3 tags:
``` python ``` python
tauADN.RTdataframe tauADN.RTdataframe
``` ```
%% Cell type:markdown id:5b80e4ef tags: %% Cell type:markdown id:5b80e4ef tags:
and the average residence time can be accessed by the attribute `RT`: and the average residence time can be accessed by the attribute `RT`:
%% Cell type:code id:c00e2b6b tags: %% Cell type:code id:c00e2b6b tags:
``` python ``` python
tauADN.RT tauADN.RT
``` ```
%% Cell type:markdown id:9dd98b49 tags: %% Cell type:markdown id:9dd98b49 tags:
The residence time can be converted in k<sub>off</sub> by: The residence time can be converted in k<sub>off</sub> by:
$k_{off}=\frac{1}{RT}$ $k_{off}=\frac{1}{RT}$
Let's create a python dictonary with the calulated k<sub>off</sub> values fro each ligand: Let's create a python dictonary with the calulated k<sub>off</sub> values fro each ligand:
%% Cell type:code id:3d0cb2d6 tags: %% Cell type:code id:3d0cb2d6 tags:
``` python ``` python
ligands = ['ADN', 'NEC'] ligands = ['ADN', 'NEC']
kinetic_parameters=[{'RL_koff':1/tauADN.RT},{'RL_koff':1/tauNEC.RT}] kinetic_parameters=[{'RL_koff':1/tauADN.RT},{'RL_koff':1/tauNEC.RT}]
``` ```
%% Cell type:markdown id:61cbf1b9 tags: %% Cell type:markdown id:61cbf1b9 tags:
## Preparation, Simulation and Analysis ## Preparation, Simulation and Analysis
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). 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).
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))* 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))*
<span style="color:black"> ⚠️ WARNING: the SSB toolkit uses μM as default concentration units.</span> <span style="color:black"> ⚠️ WARNING: the SSB toolkit uses μM as default concentration units.</span>
%% Cell type:code id:d8030846 tags: %% Cell type:code id:d8030846 tags:
``` python ``` python
# setting the range of ligand concentration # setting the range of ligand concentration
lig_conc_min = 1E-4 # μM lig_conc_min = 1E-4 # μM
lig_conc_max = 1E4 # μM lig_conc_max = 1E4 # μM
lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values lig_conc_range = np.geomspace(lig_conc_min, lig_conc_max, 20) # 20 concentration values
# Setting receptor concentration # Setting receptor concentration
receptor_conc = 1E-3 #μM receptor_conc = 1E-3 #μM
# defining other simulation parameters: # defining other simulation parameters:
time = 10000 # time of simulation in seconds time = 10000 # time of simulation in seconds
nsteps = 1000 # number of time steps nsteps = 1000 # number of time steps
``` ```
%% Cell type:markdown id:46a7f3f1 tags: %% Cell type:markdown id:46a7f3f1 tags:
## Selection of the signaling pathway ## Selection of the signaling pathway
The core of the SSB framework is, naturally, the mathematical models of the GPCRs' signaling pathways. The core of the SSB framework is, naturally, the mathematical models of the GPCRs' signaling pathways.
Since G-protein sub-families are classified by their $\alpha$ subunits, this classfication as been served to identify their signaling pathways: Since G-protein sub-families are classified by their $\alpha$ subunits, this classfication as been served to identify their signaling pathways:
* G<sub>s</sub> * G<sub>s</sub>
* G<sub>i/o</sub> * G<sub>i/o</sub>
* G<sub>q/11</sub> * G<sub>q/11</sub>
* G<sub>12/13</sub> * G<sub>12/13</sub>
📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details. 📖 See [The SSB toolkit](https://github.com/rribeiro-sci/SSBtoolkit/blob/main/docs/ssb_toolkit.md) documentation for more details.
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). 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 id:4db29d11 tags: %% Cell type:code id:4db29d11 tags:
``` python ``` python
# Selecting the signaling pathway # Selecting the signaling pathway
uniprotID = 'P29274' # Adenosine A2 receptor Uniprot ID uniprotID = 'P29274' # Adenosine A2 receptor Uniprot ID
gprotein=get.gprotein(uniprotID) gprotein=get.gprotein(uniprotID)
gprotein gprotein
``` ```
%% Cell type:markdown id:9b79bb4e tags: %% Cell type:markdown id:9b79bb4e tags:
## Integration of ODEs ## Integration of ODEs
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. 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.
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. 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.
%% Cell type:code id:5b726eb7 tags: %% Cell type:code id:5b726eb7 tags:
``` python ``` python
#create a simulation instance. #create a simulation instance.
sim = simulation.activation() sim = simulation.activation()
``` ```
%% Cell type:code id:a91aeaea tags: %% Cell type:code id:a91aeaea tags:
``` python ``` python
#Setting the simulation parameters #Setting the simulation parameters
sim.SetSimulationParameters(ligands=ligands, sim.SetSimulationParameters(ligands=ligands,
pathway=gprotein, pathway=gprotein,
receptor_conc=receptor_conc, receptor_conc=receptor_conc,
lig_conc_range=lig_conc_range, lig_conc_range=lig_conc_range,
ttotal=time, ttotal=time,
nsteps=nsteps, nsteps=nsteps,
binding_kinetics=True, binding_kinetics=True,
binding_kinetic_parameters=kinetic_parameters) binding_kinetic_parameters=kinetic_parameters)
``` ```
%% Cell type:code id:6a30d341 tags: %% Cell type:code id:6a30d341 tags:
``` python ``` python
#Running the simulation #Running the simulation
sim.Run() sim.Run()
``` ```
%% Cell type:markdown id:39d2dd3c tags: %% Cell type:markdown id:39d2dd3c tags:
In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance. In the end, the concentration values of the species of the signaling pathway over the simulation time will be saved inside the instance.
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. 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.
To achieve this, we will use the analysis attribute: To achieve this, we will use the analysis attribute:
%% Cell type:code id:d9b59f02 tags: %% Cell type:code id:d9b59f02 tags:
``` python ``` python
sim.Analysis() sim.Analysis()
``` ```
%% Cell type:markdown id:ebca39a4 tags: %% Cell type:markdown id:ebca39a4 tags:
We can now to plot the dose-response curves: We can now to plot the dose-response curves:
%% Cell type:code id:53bfafb5 tags: %% Cell type:code id:53bfafb5 tags:
``` python ``` python
sim.Curve() sim.Curve()
``` ```
%% Cell type:markdown id:f7a9ec3d tags: %% Cell type:markdown id:f7a9ec3d tags:
Finnaly, from the dose-response curves we can interpolate the EC<sub>50</sub> values. Finnaly, from the dose-response curves we can interpolate the EC<sub>50</sub> values.
%% Cell type:code id:f8300a29 tags: %% Cell type:code id:f8300a29 tags:
``` python ``` python
sim.Potency() sim.Potency()
``` ```
%% Cell type:markdown id:37ce7d4e tags: %% Cell type:markdown id:37ce7d4e tags:
💡 The potency predicted values can be exported as a python dictionary using the function `sim.PotencyToDict()` or saved in a csv file: `sim.PotencyToCSV()`. 💡 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 tags: %% Cell type:markdown id:e9e16ef7 tags:
## Save results on Google Drive ## Save results on Google Drive
%% Cell type:code id:d3ace9c6 tags: %% Cell type:code id:d3ace9c6 tags:
``` python ``` python
from google.colab import drive from google.colab import drive
drive.mount('/gdrive') drive.mount('/gdrive')
``` ```
%% Cell type:code id:2df68659 tags: %% Cell type:code id:2df68659 tags:
``` python ``` python
sim.PotencyToCSV(path='/gdrive/MyDrive/XX') ## change XX accordingly. sim.PotencyToCSV(path='/gdrive/MyDrive/XX') ## change XX accordingly.
``` ```
%% Cell type:markdown id:65febe37 tags: %% Cell type:markdown id:65febe37 tags:
## Congratulations! ## Congratulations!
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. 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.
## Cite Us ## Cite Us
If you use or adapt this tutorial for your own research projects please cite us. If you use or adapt this tutorial for your own research projects please cite us.
``` ```
@article{XXX, @article{XXX,
title={XXX}, title={XXX},
author={XXX}, author={XXX},
publisher={XXX}, publisher={XXX},
note={\url{XXX}}, note={\url{XXX}},
year={XXX} year={XXX}
} }
``` ```
## Acknowledgments ## Acknowledgments
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). 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).
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg" width="300" align='left' style="margin-left:50px">
</div> </div>
......
%% Cell type:markdown id:e532cf03 tags: %% Cell type:markdown id:e532cf03 tags:
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<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"/> <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"/>
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px">
<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"> <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">
</div> </div>
<br><br><br><br> <br><br><br><br>
<br> <br>
# Exploring SSB pathways associated to disease variants # Exploring SSB pathways associated to disease variants
## The OXTR pathway as a study case ## The OXTR pathway as a study case
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. 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.
The data used in this notebook is based on: The data used in this notebook is based on:
<blockquote> <blockquote>
Meyer, M., Jurek, B., Alfonso-Prieto, M. et al. Structure-function relationships of the Meyer, M., Jurek, B., Alfonso-Prieto, M. et al. Structure-function relationships of the
disease-linked A218T oxytocin receptor variant. Mol Psychiatry (2022). disease-linked A218T oxytocin receptor variant. Mol Psychiatry (2022).
<a href="https://doi.org/10.1038/s41380-021-01241-8" target="_bank"> https://doi.org/10.1038/s41380-021-01241-8</a> <a href="https://doi.org/10.1038/s41380-021-01241-8" target="_bank"> https://doi.org/10.1038/s41380-021-01241-8</a>
</blockquote> </blockquote>
## <span style="color:red">This tutorial just works specifically on EBRAINS</span> ## <span style="color:red">This tutorial just works specifically on EBRAINS</span>
%% Cell type:markdown id:199d4c4e tags: %% Cell type:markdown id:199d4c4e tags:
## Install SSB toolkit on EBRAINS ## Install SSB toolkit on EBRAINS
#### If you didn't clone the repository into ebrains run the following code #### If you didn't clone the repository into ebrains run the following code
%% Cell type:code id:626cb137 tags: %% Cell type:code id:626cb137 tags:
``` python ``` python
#@title Install dependencies #@title Install dependencies
%%bash %%bash
# download SSB source code and install dependencies # download SSB source code and install dependencies
if [ ! -d "SSBtoolkit/" ]; then if [ ! -d "SSBtoolkit/" ]; then
git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet git clone https://github.com/rribeiro-sci/SSBtoolkit.git --quiet
pip install -r SSBtoolkit/requirements.txt pip install -r SSBtoolkit/requirements.txt
fi fi
``` ```
%% Cell type:markdown id:42de6245 tags: %% Cell type:markdown id:42de6245 tags:
#### If you have cloned the repository run the following code #### If you have cloned the repository run the following code
%% Cell type:code id:21fadd5a tags: %% Cell type:code id:21fadd5a tags:
``` python ``` python
from Ipython.display import clear_ouput from Ipython.display import clear_output
pip install -r requirements.txt !pip install -r requirements.txt
clear_ouput() clear_ouput()
``` ```
%% Cell type:markdown id:439d9e19 tags: %% Cell type:markdown id:439d9e19 tags:
### Import Dependencies ### Import Dependencies
%% Cell type:code id:276f8a5a tags: %% Cell type:code id:276f8a5a tags:
``` python ``` python
#Import Python dependencies #Import Python dependencies
import os, sys, json, math, site import os, sys, json, math, site
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import seaborn as sns import seaborn as sns
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression from sklearn.linear_model import LinearRegression
from src.lib.ssbtoolkit import convert, get, simulation from src.lib.ssbtoolkit import convert, get, simulation
``` ```
%% Cell type:code id:745c0327 tags: %% Cell type:code id:745c0327 tags:
``` python ``` python
#Setting BioNetGen environment path #Setting BioNetGen environment path
distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/' distpath = '/opt/app-root/src/.local/lib/python3.8/site-packages/'
BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:') BioNetGen=os.path.join(distpath, 'bionetgen/bng-linux:')
mypath=%env PATH mypath=%env PATH
newpath=BioNetGen+mypath newpath=BioNetGen+mypath
%env PATH=$newpath %env PATH=$newpath
``` ```
%% Cell type:code id:970d01a1 tags: %% Cell type:code id:970d01a1 tags:
``` python ``` python
#Test bioservices #Test bioservices
from bioservices import UniProt from bioservices import UniProt
u = UniProt(verbose=False) u = UniProt(verbose=False)
``` ```
%% Cell type:markdown id:55b54990 tags: %% Cell type:markdown id:55b54990 tags:
## Introduction ## Introduction
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). 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).
The OXT exerts its function by binding to its target receptor (OXTR), a member of the Class A G-protein coupled receptors (GPCR) family. The OXT exerts its function by binding to its target receptor (OXTR), a member of the Class A G-protein coupled receptors (GPCR) family.
The mutation of an alanine residue at position 218 to a threonine on the OXTR was associated with a fenotype related to ASD. The mutation of an alanine residue at position 218 to a threonine on the OXTR was associated with a fenotype related to ASD.
In vitro experiments have shown a difference on the dynamics of the intracellular Ca<sup>2+</sup> between the OXTR-WT and OXTR-A218T. In vitro experiments have shown a difference on the dynamics of the intracellular Ca<sup>2+</sup> between the OXTR-WT and OXTR-A218T.
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. 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.
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. 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.
## Experimental Data ## Experimental Data
<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>). <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>).
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>). 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>).
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>). 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>).
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>). 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>).
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1640172296/SSBColab/OXT_results_vsvpns.png" width=600/> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1640172296/SSBColab/OXT_results_vsvpns.png" width=600/>
</div> </div>
</blockquote> </blockquote>
%% Cell type:markdown id:16404990 tags: %% Cell type:markdown id:16404990 tags:
## Modeling of the Ca<sup>2+</sup> concentration differences between WT and A218T OXTR ## Modeling of the Ca<sup>2+</sup> concentration differences between WT and A218T OXTR
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>). 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>).
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. 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.
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. 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.
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. 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 tags: %% Cell type:markdown id:8bc1a2ab tags:
At first place we need to set up the following parameters: At first place we need to set up the following parameters:
* **fitting parameters** * **fitting parameters**
* **simulation parameters** * **simulation parameters**
* **pathway parameters** * **pathway parameters**
#### Fitting Parameters #### Fitting Parameters
* expratio (flt) --> experimental signalling specie concentration ratio * expratio (flt) --> experimental signalling specie concentration ratio
* target_parameter (str) --> kinetic parameter to be modified * target_parameter (str) --> kinetic parameter to be modified
* maxiter (int) --> maximum number of iteration * maxiter (int) --> maximum number of iteration
* seed (flt) --> ramdom seed for scaling the modified parameter * seed (flt) --> ramdom seed for scaling the modified parameter
* seed_incrementor (flt) --> seed incrementor (each iteration will increment the seed by this value) * seed_incrementor (flt) --> seed incrementor (each iteration will increment the seed by this value)
* seed_decrementor (flt) --> seed decrementor (each iteration will decrement the seed by this value) * seed_decrementor (flt) --> seed decrementor (each iteration will decrement the seed by this value)
#### Simulation Parameters #### Simulation Parameters
* ttotal (int) --> time of simulation in seconds * ttotal (int) --> time of simulation in seconds
* nsteps (int) --> number of time steps * nsteps (int) --> number of time steps
* pathway (str) --> name of the signalling pathway * pathway (str) --> name of the signalling pathway
* observable (str) --> signalling specie to be measured * observable (str) --> signalling specie to be measured
#### Pahtway Parameters #### Pahtway Parameters
* time_in (int) --> simulation time at wich the equation will start to be simulated * time_in (int) --> simulation time at wich the equation will start to be simulated
* time_out (int) --> simulation time at wich the equation will stop to be simulated * time_out (int) --> simulation time at wich the equation will stop to be simulated
%% Cell type:code id:b007dfa0 tags: %% Cell type:code id:b007dfa0 tags:
``` python ``` python
# we will start by creating a simulation instance. # we will start by creating a simulation instance.
sim = simulation.fitModel() sim = simulation.fitModel()
``` ```
%% Cell type:code id:4e1d1b92 tags: %% Cell type:code id:4e1d1b92 tags:
``` python ``` python
# configuration of simulation parameters # configuration of simulation parameters
sim.SetSimulationParameters(ttotal=250, sim.SetSimulationParameters(ttotal=250,
nsteps=10000, nsteps=10000,
pathway='OXTR_pathway', pathway='OXTR_pathway',
observable='Ca', observable='Ca',
pathway_parameters={'time_in':50, pathway_parameters={'time_in':50,
'time_out':51}) 'time_out':51})
``` ```
%% Cell type:code id:4cc4a960 tags: %% Cell type:code id:4cc4a960 tags:
``` python ``` python
#Start the simulation #Start the simulation
sim.Run(expratio=1.13, sim.Run(expratio=1.13,
target_parameter='R_L_Gq_trimer_kf', target_parameter='R_L_Gq_trimer_kf',
maxiter=20, maxiter=20,
seed=0.2, seed=0.2,
seed_incrementor=0.5, seed_incrementor=0.5,
seed_decrementor=0.1) seed_decrementor=0.1)
``` ```
%% Cell type:markdown id:bd653d12 tags: %% Cell type:markdown id:bd653d12 tags:
We can now plot the curves, amplitudes, area under the curve (AUC), and the full width at half maximum (FWHM): We can now plot the curves, amplitudes, area under the curve (AUC), and the full width at half maximum (FWHM):
%% Cell type:code id:d37084ed tags: %% Cell type:code id:d37084ed tags:
``` python ``` python
sim.plotCurves() sim.plotCurves()
``` ```
%% Cell type:markdown id:5fd46fe3 tags: %% Cell type:markdown id:5fd46fe3 tags:
## Save results on Google Drive ## Save results on Google Drive
%% Cell type:code id:3e476cad tags: %% Cell type:code id:3e476cad tags:
``` python ``` python
from google.colab import drive from google.colab import drive
drive.mount('/gdrive') drive.mount('/gdrive')
``` ```
%% Cell type:code id:b98afd3b tags: %% Cell type:code id:b98afd3b tags:
``` python ``` python
sim.plotCurves(save=True, filename='/gdrive/MyDrive/XX') ## change XX accordingly. sim.plotCurves(save=True, filename='/gdrive/MyDrive/XX') ## change XX accordingly.
``` ```
%% Cell type:markdown id:ff3d7f41 tags: %% Cell type:markdown id:ff3d7f41 tags:
## Congratulations! ## Congratulations!
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. 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.
## Cite Us ## Cite Us
If you use or adapt this tutorial for your own research projects please cite us. If you use or adapt this tutorial for your own research projects please cite us.
``` ```
@article{XXX, @article{XXX,
title={XXX}, title={XXX},
author={XXX}, author={XXX},
publisher={XXX}, publisher={XXX},
note={\url{XXX}}, note={\url{XXX}},
year={XXX} year={XXX}
} }
``` ```
## Acknowledgments ## Acknowledgments
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). 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).
<div style="padding-bottom:50px"> <div style="padding-bottom:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1637657234/logos/HBP_horizontal_logo_qtcyzn.png" width="300" align='left' style="margin-left:50px">
<img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg" width="300" align='left' style="margin-left:50px"> <img src="https://res.cloudinary.com/djz27k5hg/image/upload/v1642677502/logos/COFUNDED_EU_j2ktlp.jpg" width="300" align='left' style="margin-left:50px">
</div> </div>
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment