|
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
class DoseResponseExperiment(SimulationExperiment):
"""Hormone dose-response curves."""
@timeit
def models(self) -> Dict[str, Union[AbstractModel, Path]]:
return {"model1": Path(__file__).parent.parent / "model" / "liver_glucose.xml"}
@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, ureg=self.ureg, udict=udict
)
return dsets
def datagenerators(self) -> None:
for key in ["time", "glu", "ins", "epi", "gamma"]:
Data(experiment=self, task="task_glc_scan", index=key)
@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=30), "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=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}