Experiments
DoseResponseExperiment
Model
Datasets
Simulations
Scans
Figures
fig1
Code
experiments/dose_response.py
from typing import Dict
import numpy as np
import pandas as pd
import xarray as xr
from sbmlsim.experiment import SimulationExperiment
from sbmlsim.data import DataSet, load_pkdb_dataframe
from sbmlsim.model import AbstractModel, RoadrunnerSBMLModel
from sbmlsim.simulation import Timecourse, TimecourseSim, ScanSim, Dimension
from sbmlsim.task import Task
from sbmlsim.utils import timeit
from matplotlib.pyplot import Figure
from sbmlsim.plot.plotting_matplotlib import add_data, plt
class DoseResponseExperiment(SimulationExperiment):
"""Hormone dose-response curves."""
@timeit
def models(self) -> Dict[str, AbstractModel]:
return {
"model1": RoadrunnerSBMLModel(source="model/liver_glucose.xml",
base_path=self.base_path)
}
@timeit
def datasets(self) -> Dict[str, DataSet]:
dsets = {}
# dose-response data for hormones
for hormone_key in ['Epinephrine', 'Glucagon', 'Insulin']:
df = load_pkdb_dataframe(
f"DoseResponse_Tab{hormone_key}",
data_path=self.data_path
)
df = df[df.condition == "normal"] # only healthy controls
epi_normal_studies = [
"Degn2004", "Lerche2009", "Mitrakou1991",
"Levy1998", "Israelian2006", "Jones1998",
"Segel2002"
]
glu_normal_studies = [
"Butler1991", "Cobelli2010", "Fery1993"
"Gerich1993", "Henkel2005",
"Mitrakou1991"
"Basu2009", "Mitrakou1992",
"Degn2004", "Lerche2009",
"Levy1998", "Israelian2006",
"Segel2002"
]
ins_normal_studies = [
'Ferrannini1988', 'Fery1993', 'Gerich1993', 'Basu2009',
'Lerche2009', 'Henkel2005', 'Butler1991', 'Knop2007',
'Cobelli2010', 'Mitrakou1992',
]
# filter studies
if hormone_key == "Epinephrine":
df = df[df.reference.isin(epi_normal_studies)]
elif hormone_key == "Glucagon":
df = df[df.reference.isin(glu_normal_studies)]
# correct glucagon data for insulin suppression
# (hyperinsulinemic clamps)
insulin_supression = 3.4
glu_clamp_studies = [
"Degn2004", "Lerche2009", "Levy1998", "Israelian2006",
"Segel2002"
]
df.loc[df.reference.isin(
glu_clamp_studies), 'mean'] = insulin_supression * df[
df.reference.isin(glu_clamp_studies)]['mean']
df.loc[df.reference.isin(
glu_clamp_studies), 'se'] = insulin_supression * df[
df.reference.isin(glu_clamp_studies)]['se']
elif hormone_key == "Insulin":
df = df[df.reference.isin(ins_normal_studies)]
udict = {
'glc': df["glc_unit"].unique()[0],
'mean': df["unit"].unique()[0],
}
dsets[hormone_key.lower()] = DataSet.from_df(df, udict=udict,
ureg=self.ureg)
return dsets
@timeit
def tasks(self) -> Dict[str, Task]:
"""Tasks"""
return {
"task_glc_scan": Task(model="model1", simulation="glc_scan")
}
@timeit
def simulations(self) -> Dict[str, ScanSim]:
"""Scanning dose-response curves of hormones and gamma function.
Vary external glucose concentrations (boundary condition).
"""
glc_scan = ScanSim(
simulation=TimecourseSim([
Timecourse(start=0, end=1, steps=1, changes={})
]),
dimensions=[
Dimension(
"dim1", changes={
'[glc_ext]': self.Q_(np.linspace(2, 20, num=10), 'mM')
}),
]
)
return {
"glc_scan": glc_scan
}
@timeit
def figures(self) -> Dict[str, Figure]:
xunit = "mM"
yunit_hormone = "pmol/l"
yunit_gamma = "dimensionless"
f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(10, 10))
f.subplots_adjust(wspace=.3, hspace=.3)
axes = (ax1, ax2, ax3, ax4)
# process scan results
task = self._tasks["task_glc_scan"]
model = self._models[task.model_id]
tcscan = self._simulations[task.simulation_id] # TimecourseScan Definition
# FIXME: this must be simpler
glc_vec = tcscan.dimensions[0].changes['[glc_ext]']
xres = self.results["task_glc_scan"] # type: Result
# we already have all the data ordered, we only want the steady state value
dose_response = {}
for sid in ["glu", "epi", "ins", "gamma"]:
da = xres[sid] # type: xr.DataArray
# get initial time
head = da.head({'_time': 1}).to_series()
dose_response[sid] = head.values
dose_response['[glc_ext]'] = glc_vec
df = pd.DataFrame(dose_response)
dset = DataSet.from_df(df, udict=model.udict, ureg=self.ureg)
print(dset)
# plot scan results
kwargs = {
'linewidth': 2,
'linestyle': '-',
'marker': 'None',
'color': 'black'
}
add_data(ax1, dset, xid="[glc_ext]", yid="glu",
xunit=xunit, yunit=yunit_hormone, **kwargs)
add_data(ax2, dset, xid="[glc_ext]", yid="epi",
xunit=xunit, yunit=yunit_hormone, **kwargs)
add_data(ax3, dset, xid="[glc_ext]", yid="ins",
xunit=xunit, yunit=yunit_hormone, **kwargs)
add_data(ax4, dset, xid="[glc_ext]", yid="gamma",
xunit=xunit, yunit=yunit_gamma, **kwargs)
# plot experimental data
kwargs = {
'color': "black",
'linestyle': "None",
'alpha': 0.6,
}
add_data(ax1, self._datasets["glucagon"], xid="glc", yid="mean",
yid_se="mean_se",
xunit=xunit, yunit=yunit_hormone, label="Glucagon", **kwargs)
add_data(ax2, self._datasets["epinephrine"], xid="glc", yid="mean",
yid_se="mean_se",
xunit=xunit, yunit=yunit_hormone, label="Epinephrine",
**kwargs)
add_data(ax3, self._datasets["insulin"], xid="glc", yid="mean",
yid_se="mean_se",
xunit=xunit, yunit=yunit_hormone, label="Insulin", **kwargs)
ax1.set_ylabel(f'glucagon [{yunit_hormone}]')
ax1.set_ylim(0, 200)
ax2.set_ylabel(f'epinephrine [{yunit_hormone}]')
ax2.set_ylim(0, 7000)
ax3.set_ylabel(f'insulin [{yunit_hormone}]')
ax3.set_ylim(0, 800)
ax4.set_ylabel(f'gamma [{yunit_gamma}]')
ax4.set_ylim(0, 1)
for ax in axes:
ax.set_xlabel(f'glucose [{xunit}]')
ax.set_xlim(2, 20)
ax2.set_xlim(2, 8)
return {
'fig1': f
}
