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_deprecated_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
@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}
[docs] def data(self) -> Dict[str, Data]: self.add_selections_data( selections=["time", "glu", "ins", "epi", "gamma"], task_ids=["task_glc_scan"], ) return {}
[docs] def figures_mpl(self) -> Dict[str, Figure]: xunit = "mM" yunit_hormone = "pmol/l" yunit_gamma = "dimensionless" fig_mpl, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(10, 10)) fig_mpl.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] # FIXME: this must be simpler glc_vec = tcscan.dimensions[0].changes["[glc_ext]"] xres: XResult = self.results["task_glc_scan"] # we already have all the data ordered, we only want the steady state value dose_response = {} for sid in ["glu", "epi", "ins", "gamma"]: da: xr.DataArray = xres[sid] # 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": fig_mpl}