Experiments

DoseResponseExperiment

Models

Datasets

Figures

fig1

Code

../../experiments/dose_response.py

from pathlib import Path
from typing import Dict

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.simulation import Dimension, ScanSim, Timecourse, TimecourseSim
from sbmlsim.task import Task
from sbmlsim.utils import timeit


class DoseResponseExperiment(SimulationExperiment):
    """Hormone dose-response curves."""

    @timeit
    def models(self) -> Dict[str, AbstractModel]:
        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, udict=udict, ureg=self.ureg
            )

        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: 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}