Experiments

DoseResponseExperiment

Model

Datasets

Simulations

Scans

Figures

DoseResponseExperiment

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
        }