Source code for sbmlsim.examples.experiments.glucose.experiments.dose_response
from pathlib import Path
from typing import Dict, Union
import numpy as np
import pandas as pd
import xarray as xr
from matplotlib.pyplot import Figure
from sbmlsim.data import Data, DataSet, load_pkdb_dataframe
from sbmlsim.experiment import SimulationExperiment
from sbmlsim.model import AbstractModel, RoadrunnerSBMLModel
from sbmlsim.plot.plotting_matplotlib import add_data, plt
from sbmlsim.result import XResult
from sbmlsim.simulation import Dimension, ScanSim, Timecourse, TimecourseSim
from sbmlsim.task import Task
from sbmlsim.units import UnitsInformation
from sbmlsim.utils import timeit
[docs]class DoseResponseExperiment(SimulationExperiment):
"""Hormone dose-response curves."""
@timeit
[docs] def models(self) -> Dict[str, Union[AbstractModel, Path]]:
return {"model1": Path(__file__).parent.parent / "model" / "liver_glucose.xml"}
@timeit
[docs] 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, ureg=self.ureg, udict=udict
)
return dsets
[docs] def datagenerators(self) -> None:
for key in ["time", "glu", "ins", "epi", "gamma"]:
Data(experiment=self, task="task_glc_scan", index=key)
@timeit
[docs] def tasks(self) -> Dict[str, Task]:
"""Tasks"""
return {"task_glc_scan": Task(model="model1", simulation="glc_scan")}
@timeit
[docs] 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=30), "mM")},
),
],
)
return {"glc_scan": glc_scan}
@timeit
[docs] 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=0.3, hspace=0.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: XResult
# 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.uinfo.udict, ureg=self.ureg)
# 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}