docs/tutorials/constructing_models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Constructing Models\n",
    "\n",
    "In this notebook example, a step-by-step approach of building a simple model$^1$ of trafficking of high-energy phosphate bonds is demonstrated. Illustrated below is a graphical view of the full system along with the reaction rate equations and numerical values:\n",
    "\n",
    "![PhosTrafficking](../media/images/phosphate_trafficking.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The example starts by creating a model of the \"use\", \"distr\", and \"form\" reactions. Then the model is expanded to include the remaining metabolites, reactions, and any additional information that should be defined in the model. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:03.439679Z",
     "iopub.status.busy": "2021-01-25T18:00:03.437958Z",
     "iopub.status.idle": "2021-01-25T18:00:05.780946Z",
     "shell.execute_reply": "2021-01-25T18:00:05.781418Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from mass import (\n",
    "    MassConfiguration, MassMetabolite, MassModel, MassReaction)\n",
    "from mass.util.matrix import left_nullspace, nullspace\n",
    "\n",
    "mass_config = MassConfiguration()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step to creating the model is to define the `MassModel` object. A `MassModel` only requires an identifier to be initialized. For best practice, it is recommended to utilize SBML compliant identifiers for all objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.791270Z",
     "iopub.status.busy": "2021-01-25T18:00:05.790612Z",
     "iopub.status.idle": "2021-01-25T18:00:05.792714Z",
     "shell.execute_reply": "2021-01-25T18:00:05.793165Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter Username\n"
     ]
    }
   ],
   "source": [
    "model = MassModel(\"Phosphate_Trafficking\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model is initially empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.797764Z",
     "iopub.status.busy": "2021-01-25T18:00:05.797099Z",
     "iopub.status.idle": "2021-01-25T18:00:05.800637Z",
     "shell.execute_reply": "2021-01-25T18:00:05.801102Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of metabolites: 0\n",
      "Number of initial conditions: 0\n",
      "Number of reactions: 0\n"
     ]
    }
   ],
   "source": [
    "print(\"Number of metabolites: {0}\".format(len(model.metabolites)))\n",
    "print(\"Number of initial conditions: {0}\".format(len(model.initial_conditions)))\n",
    "print(\"Number of reactions: {0}\".format(len(model.reactions)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to create `MassMetabolite` and `MassReaction` objects to represent the metabolites and reactions that should exist in the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining metabolites\n",
    "\n",
    "To create a `MassMetabolite`, a unique identifier is required. The `formula` and `charge` attributes are set to ensure mass and charge balancing of reactions in which the metabolite is a participant. The `compartment` attribute indicates where the metabolite is located. In this model, all metabolites exist in a single compartment, abbreviated as \"c\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.806713Z",
     "iopub.status.busy": "2021-01-25T18:00:05.805996Z",
     "iopub.status.idle": "2021-01-25T18:00:05.808067Z",
     "shell.execute_reply": "2021-01-25T18:00:05.808512Z"
    }
   },
   "outputs": [],
   "source": [
    "atp_c = MassMetabolite(\n",
    "    \"atp_c\",\n",
    "    name=\"ATP\",\n",
    "    formula=\"C10H12N5O13P3\",\n",
    "    charge=-4,\n",
    "    compartment=\"c\")\n",
    "\n",
    "adp_c = MassMetabolite(\n",
    "    \"adp_c\",\n",
    "    name=\"ADP\",\n",
    "    formula=\"C10H12N5O10P2\",\n",
    "    charge=-3,\n",
    "    compartment=\"c\")\n",
    "\n",
    "amp_c = MassMetabolite(\n",
    "    \"amp_c\",\n",
    "    name=\"AMP\",\n",
    "    formula=\"C10H12N5O7P\",\n",
    "    charge=-2,\n",
    "    compartment=\"c\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The metabolite concentrations can be defined as the initial conditions for the metabolites using the `initial_condition` attribute. As previously stated, the concentrations are $\\text{[ATP]}=1.6$, $\\text{[ADP]}=0.4$, and $\\text{[AMP]}=0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.813223Z",
     "iopub.status.busy": "2021-01-25T18:00:05.812519Z",
     "iopub.status.idle": "2021-01-25T18:00:05.815239Z",
     "shell.execute_reply": "2021-01-25T18:00:05.815693Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "atp_c: 1.6\n",
      "adp_c: 0.4\n",
      "amp_c: 0.1\n"
     ]
    }
   ],
   "source": [
    "atp_c.initial_condition = 1.6\n",
    "adp_c.ic = 0.4  # Alias for initial_condition\n",
    "amp_c.ic = 0.1\n",
    "\n",
    "for metabolite in [atp_c, adp_c, amp_c]:\n",
    "    print(\"{0}: {1}\".format(metabolite.id, metabolite.initial_condition))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The metabolites are currently not a part of any reaction. Consequently, the `ordinary_differential_equation` attribute is `None`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.820604Z",
     "iopub.status.busy": "2021-01-25T18:00:05.819478Z",
     "iopub.status.idle": "2021-01-25T18:00:05.822590Z",
     "shell.execute_reply": "2021-01-25T18:00:05.823038Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "None\n",
      "None\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "print(atp_c.ordinary_differential_equation)\n",
    "print(adp_c.ode)  # Alias for ordinary_differential_equation\n",
    "print(amp_c.ode)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to create the reactions in which the metabolites participate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining reactions\n",
    "\n",
    "Just like `MassMetabolite` objects, a unique identifier is also required to create a `MassReaction`. The `reversible` attribute determines whether the reaction can proceed in both the forward and reverse directions, or only in the forward direction. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.827789Z",
     "iopub.status.busy": "2021-01-25T18:00:05.827100Z",
     "iopub.status.idle": "2021-01-25T18:00:05.829148Z",
     "shell.execute_reply": "2021-01-25T18:00:05.829721Z"
    }
   },
   "outputs": [],
   "source": [
    "distr = MassReaction(\"distr\", name=\"Distribution\", reversible=True)\n",
    "use = MassReaction(\"use\", name=\"ATP Utilization\", reversible=False)\n",
    "form = MassReaction(\"form\", name=\"ATP Formation\", reversible=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metabolites are added to reactions using a dictionary of metabolite objects and their stoichiometric coefficients. A group of metabolites can be added either all at once or one at a time. A negative coefficient indicates the metabolite is a reactant, while a positive coefficient indicates the metabolite is a product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.835157Z",
     "iopub.status.busy": "2021-01-25T18:00:05.834411Z",
     "iopub.status.idle": "2021-01-25T18:00:05.837220Z",
     "shell.execute_reply": "2021-01-25T18:00:05.837713Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "distr: 2 adp_c <=> amp_c + atp_c\n",
      "use: atp_c --> adp_c\n",
      "form: adp_c --> atp_c\n"
     ]
    }
   ],
   "source": [
    "distr.add_metabolites({\n",
    "    adp_c: -2,\n",
    "    atp_c: 1,\n",
    "    amp_c: 1})\n",
    "\n",
    "use.add_metabolites({\n",
    "    atp_c: -1,\n",
    "    adp_c: 1})\n",
    "\n",
    "form.add_metabolites({\n",
    "    adp_c: -1,\n",
    "    atp_c: 1})\n",
    "\n",
    "for reaction in [distr, use, form]:\n",
    "    print(reaction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the reactions are created, their parameters can be defined. As stated earlier, the distribution reaction is considerably faster when compared to other reactions in the model. The forward rate constant $k^{\\rightarrow}$, represented as `kf`, can be set as $k^{\\rightarrow}_{distr}=1000\\ \\text{min}^{-1}$. The equilibrium constant $K_{eq}$, represented as `Keq`, is approximately $K_{distr}=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.844635Z",
     "iopub.status.busy": "2021-01-25T18:00:05.843887Z",
     "iopub.status.idle": "2021-01-25T18:00:05.846838Z",
     "shell.execute_reply": "2021-01-25T18:00:05.847287Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'kf_distr': 1000, 'Keq_distr': 1}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distr.forward_rate_constant = 1000\n",
    "distr.equilibrium_constant = 1\n",
    "distr.parameters  # Return defined mass action kinetic parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown earlier, the forward rate constants are set as $k^{\\rightarrow}_{use}=6.25\\ \\text{min}^{-1}$ and $k^{\\rightarrow}_{form}=25\\ \\text{min}^{-1}$.\n",
    "The `kf_str` attribute can be used to get the identifier of the forward rate constant as a string."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.852526Z",
     "iopub.status.busy": "2021-01-25T18:00:05.851802Z",
     "iopub.status.idle": "2021-01-25T18:00:05.854494Z",
     "shell.execute_reply": "2021-01-25T18:00:05.855015Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kf_use: 6.25\n",
      "kf_form: 25\n"
     ]
    }
   ],
   "source": [
    "use.forward_rate_constant = 6.25\n",
    "form.kf = 25  # Alias for forward_rate_constant\n",
    "\n",
    "print(\"{0}: {1}\".format(use.kf_str, use.kf))\n",
    "print(\"{0}: {1}\".format(form.kf_str, form.kf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reactions can be added to the model using the `add_reactions()` method. Adding the reactions to the model also adds the associated metabolites and genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.862410Z",
     "iopub.status.busy": "2021-01-25T18:00:05.861670Z",
     "iopub.status.idle": "2021-01-25T18:00:05.864439Z",
     "shell.execute_reply": "2021-01-25T18:00:05.864894Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of metabolites: 3\n",
      "Number of initial conditions: 3\n",
      "Number of reactions: 3\n"
     ]
    }
   ],
   "source": [
    "model.add_reactions([distr, use, form])\n",
    "\n",
    "print(\"Number of metabolites: {0}\".format(len(model.metabolites)))\n",
    "print(\"Number of initial conditions: {0}\".format(len(model.initial_conditions)))\n",
    "print(\"Number of reactions: {0}\".format(len(model.reactions)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stoichiometric matrix of the model is automatically constructed with the addition of the reactions and metabolites to the model. It can be accessed through the `stoichiometric_matrix` property (alias `S`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.869754Z",
     "iopub.status.busy": "2021-01-25T18:00:05.868951Z",
     "iopub.status.idle": "2021-01-25T18:00:05.872077Z",
     "shell.execute_reply": "2021-01-25T18:00:05.872709Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-2.  1. -1.]\n",
      " [ 1. -1.  1.]\n",
      " [ 1.  0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "print(model.S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stoichiometric matrix attribute can be updated and stored in various formats using the `update_S()` method. For example, the stoichiometric matrix can be converted and stored as a `pandas.DataFrame`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.879779Z",
     "iopub.status.busy": "2021-01-25T18:00:05.878935Z",
     "iopub.status.idle": "2021-01-25T18:00:05.887582Z",
     "shell.execute_reply": "2021-01-25T18:00:05.888145Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>distr</th>\n",
       "      <th>use</th>\n",
       "      <th>form</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>adp_c</th>\n",
       "      <td>-2</td>\n",
       "      <td>1</td>\n",
       "      <td>-1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>atp_c</th>\n",
       "      <td>1</td>\n",
       "      <td>-1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>amp_c</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       distr  use  form\n",
       "adp_c     -2    1    -1\n",
       "atp_c      1   -1     1\n",
       "amp_c      1    0     0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.update_S(array_type=\"DataFrame\", dtype=np.int_, update_model=True)\n",
    "model.S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Associating the metabolites with reactions allows for the mass action reaction rate expressions to be generated based on the stoichiometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.908047Z",
     "iopub.status.busy": "2021-01-25T18:00:05.907106Z",
     "iopub.status.idle": "2021-01-25T18:00:05.910200Z",
     "shell.execute_reply": "2021-01-25T18:00:05.910700Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kf_distr*(adp_c(t)**2 - amp_c(t)*atp_c(t)/Keq_distr)\n"
     ]
    }
   ],
   "source": [
    "print(distr.rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generation of the reaction rates also allows for the metabolite ODEs to be generated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.925050Z",
     "iopub.status.busy": "2021-01-25T18:00:05.924275Z",
     "iopub.status.idle": "2021-01-25T18:00:05.927338Z",
     "shell.execute_reply": "2021-01-25T18:00:05.927861Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kf_distr*(adp_c(t)**2 - amp_c(t)*atp_c(t)/Keq_distr) + kf_form*adp_c(t) - kf_use*atp_c(t)\n"
     ]
    }
   ],
   "source": [
    "print(atp_c.ode)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `nullspace()` method can be used to obtain the null space of the stoichiometric matrix. The nullspace reflects the pathways through the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.937717Z",
     "iopub.status.busy": "2021-01-25T18:00:05.936544Z",
     "iopub.status.idle": "2021-01-25T18:00:05.940063Z",
     "shell.execute_reply": "2021-01-25T18:00:05.940563Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Pathway 1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>distr</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>use</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>form</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Pathway 1\n",
       "distr          0\n",
       "use            1\n",
       "form           1"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ns = nullspace(model.S)  # Get the null space\n",
    "# Divide by the minimum and round to nearest integer\n",
    "ns = np.rint(ns / np.min(ns[np.nonzero(ns)]))\n",
    "pd.DataFrame(\n",
    "    ns, index=model.reactions.list_attr(\"id\"),  # Rows represent reactions\n",
    "    columns=[\"Pathway 1\"], dtype=np.int_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a similar fashion the left nullspace can be obtained using the `left_nullspace` function. The left nullspace represents the conserved moieties in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.951692Z",
     "iopub.status.busy": "2021-01-25T18:00:05.950293Z",
     "iopub.status.idle": "2021-01-25T18:00:05.957619Z",
     "shell.execute_reply": "2021-01-25T18:00:05.958110Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>adp_c</th>\n",
       "      <th>atp_c</th>\n",
       "      <th>amp_c</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Total AxP</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           adp_c  atp_c  amp_c\n",
       "Total AxP      1      1      1"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lns = left_nullspace(model.S)\n",
    "# Divide by the minimum and round to nearest integer\n",
    "lns = np.rint(lns / np.min(lns[np.nonzero(lns)]))\n",
    "pd.DataFrame(\n",
    "    lns, index=[\"Total AxP\"],\n",
    "    columns=model.metabolites.list_attr(\"id\"),  # Columns represent metabolites\n",
    "    dtype=np.int_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expanding an Existing Model\n",
    "\n",
    "Now, the existing model is expanded to include a buffer reaction, where a phosphagen is utilized to store a high-energy phosphate in order to buffer the ATP/ADP ratio as needed. Because the buffer molecule represents a generic phosphagen, there is no chemical formula for the molecule. Therefore, the buffer molecule can be represented as a moiety in the `formula` attribute using square brackets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.967852Z",
     "iopub.status.busy": "2021-01-25T18:00:05.966583Z",
     "iopub.status.idle": "2021-01-25T18:00:05.970036Z",
     "shell.execute_reply": "2021-01-25T18:00:05.970952Z"
    }
   },
   "outputs": [],
   "source": [
    "b = MassMetabolite(\n",
    "    \"B\",\n",
    "    name=\"Phosphagen buffer (Free)\",\n",
    "    formula=\"[B]\",\n",
    "    charge=0,\n",
    "    compartment=\"c\")\n",
    "\n",
    "bp = MassMetabolite(\n",
    "    \"BP\",\n",
    "    name=\"Phosphagen buffer (Loaded)\",\n",
    "    formula=\"[B]-PO3\",\n",
    "    charge=-1,\n",
    "    compartment=\"c\")\n",
    "\n",
    "buffer = MassReaction(\"buffer\", name=\"ATP Buffering\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When adding metabolites to the reaction, the `get_by_id()` method is used to add already existing metabolites in the model to the reaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.984417Z",
     "iopub.status.busy": "2021-01-25T18:00:05.983027Z",
     "iopub.status.idle": "2021-01-25T18:00:05.986065Z",
     "shell.execute_reply": "2021-01-25T18:00:05.986955Z"
    }
   },
   "outputs": [],
   "source": [
    "buffer.add_metabolites({\n",
    "    b: -1,\n",
    "    model.metabolites.get_by_id(\"atp_c\"): -1,\n",
    "    model.metabolites.get_by_id(\"adp_c\"): 1,\n",
    "    bp: 1})\n",
    "\n",
    "# Add reaction to model\n",
    "model.add_reactions([buffer])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this reaction, $k^{\\rightarrow}_{buffer}=1000\\ \\text{min}^{-1}$ and $K_{buffer}=1$. Because the reaction has already been added to the model, the `MassModel.update_parameters()` method can be used to update the reaction parameters using a dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:05.996046Z",
     "iopub.status.busy": "2021-01-25T18:00:05.994787Z",
     "iopub.status.idle": "2021-01-25T18:00:05.999536Z",
     "shell.execute_reply": "2021-01-25T18:00:06.000616Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'kf_buffer': 1000, 'Keq_buffer': 1}"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.update_parameters({\n",
    "    buffer.kf_str: 1000,\n",
    "    buffer.Keq_str: 1})\n",
    "\n",
    "buffer.parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By adding the reaction to the model, the left nullspace expanded to include a conservation pool for the total buffer in the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.014397Z",
     "iopub.status.busy": "2021-01-25T18:00:06.013543Z",
     "iopub.status.idle": "2021-01-25T18:00:06.016519Z",
     "shell.execute_reply": "2021-01-25T18:00:06.016953Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>adp_c</th>\n",
       "      <th>atp_c</th>\n",
       "      <th>amp_c</th>\n",
       "      <th>B</th>\n",
       "      <th>BP</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Total AxP</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Total Buffer</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              adp_c  atp_c  amp_c  B  BP\n",
       "Total AxP         1      1      1  0   0\n",
       "Total Buffer      0      0      0  1   1"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lns = left_nullspace(model.S)\n",
    "for i, row in enumerate(lns):\n",
    "    # Divide by the minimum and round to nearest integer\n",
    "    lns[i] = np.rint(row / np.min(row[np.nonzero(row)]))\n",
    "pd.DataFrame(lns, index=[\"Total AxP\", \"Total Buffer\"],\n",
    "             columns=model.metabolites.list_attr(\"id\"),\n",
    "             dtype=np.int_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performing symbolic calculations\n",
    "\n",
    "Although the concentrations for the free and loaded buffer molecules are currently unknown, the total amount of buffer is known and set as $B_{total} = 10$. Because the buffer reaction is assumed to be at equilibrium, it becomes possible to solve for the concentrations of the free and loaded buffer molecules.\n",
    "\n",
    "Below, the symbolic capabilities of **SymPy** are used to solve for the steady state concentrations of the buffer molecules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.021706Z",
     "iopub.status.busy": "2021-01-25T18:00:06.020808Z",
     "iopub.status.idle": "2021-01-25T18:00:06.023583Z",
     "shell.execute_reply": "2021-01-25T18:00:06.023052Z"
    }
   },
   "outputs": [],
   "source": [
    "from sympy import Eq, Symbol, pprint, solve\n",
    "\n",
    "from mass.util import strip_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to define the equation for the total buffer pool symbolically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.036470Z",
     "iopub.status.busy": "2021-01-25T18:00:06.035738Z",
     "iopub.status.idle": "2021-01-25T18:00:06.038710Z",
     "shell.execute_reply": "2021-01-25T18:00:06.039397Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "B + BP = 10\n"
     ]
    }
   ],
   "source": [
    "buffer_total_equation = Eq(Symbol(\"B\") + Symbol(\"BP\"), 10)\n",
    "pprint(buffer_total_equation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation for the reaction rate at equilibrium is also defined. The `strip_time()` function is used to strip time dependency from the equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.090869Z",
     "iopub.status.busy": "2021-01-25T18:00:06.088650Z",
     "iopub.status.idle": "2021-01-25T18:00:06.093172Z",
     "shell.execute_reply": "2021-01-25T18:00:06.093682Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 = 1600.0⋅B - 400.0⋅BP\n"
     ]
    }
   ],
   "source": [
    "buffer_rate_equation = Eq(0, strip_time(buffer.rate))\n",
    "# Substitute defined concentration values into equation\n",
    "buffer_rate_equation = buffer_rate_equation.subs({\n",
    "    \"atp_c\":  atp_c.initial_condition,\n",
    "    \"adp_c\":  adp_c.initial_condition,\n",
    "    \"kf_buffer\": buffer.kf,\n",
    "    \"Keq_buffer\": buffer.Keq})\n",
    "pprint(buffer_rate_equation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These two equations can be solved to get the buffer concentrations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.116101Z",
     "iopub.status.busy": "2021-01-25T18:00:06.115252Z",
     "iopub.status.idle": "2021-01-25T18:00:06.118110Z",
     "shell.execute_reply": "2021-01-25T18:00:06.118561Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{B: 2.00000000000000, BP: 8.00000000000000}"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "buffer_sol = solve([buffer_total_equation, buffer_rate_equation],\n",
    "                   [Symbol(\"B\"), Symbol(\"BP\")])\n",
    "buffer_sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the metabolites already exist in the model, their initial conditions can be updated to the calculated concentrations using the `MassModel.update_initial_conditions()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.124344Z",
     "iopub.status.busy": "2021-01-25T18:00:06.123558Z",
     "iopub.status.idle": "2021-01-25T18:00:06.126475Z",
     "shell.execute_reply": "2021-01-25T18:00:06.127068Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{<MassMetabolite adp_c at 0x7f8368c373a0>: 0.4,\n",
       " <MassMetabolite atp_c at 0x7f8368c377c0>: 1.6,\n",
       " <MassMetabolite amp_c at 0x7f8368c37f70>: 0.1,\n",
       " <MassMetabolite B at 0x7f8359672160>: 2.0,\n",
       " <MassMetabolite BP at 0x7f8359672b50>: 8.0}"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Replace the symbols in the dict\n",
    "for met_symbol, concentration in buffer_sol.copy().items():\n",
    "    metabolite = model.metabolites.get_by_id(str(met_symbol))\n",
    "    # Make value as a float\n",
    "    buffer_sol[metabolite] = float(buffer_sol.pop(met_symbol))\n",
    "\n",
    "model.update_initial_conditions(buffer_sol)\n",
    "model.initial_conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding boundary reactions\n",
    "\n",
    "After adding the buffer reactions, the next step is to define the AMP source and demand reactions. The `add_boundary()` method is employed to create and add a boundary reaction to a model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.133880Z",
     "iopub.status.busy": "2021-01-25T18:00:06.133164Z",
     "iopub.status.idle": "2021-01-25T18:00:06.135796Z",
     "shell.execute_reply": "2021-01-25T18:00:06.136292Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "amp_drain: amp_c --> \n",
      "amp_in: amp_c --> \n"
     ]
    }
   ],
   "source": [
    "amp_drain = model.add_boundary(\n",
    "    model.metabolites.amp_c,\n",
    "    boundary_type=\"demand\",\n",
    "    reaction_id=\"amp_drain\")\n",
    "\n",
    "amp_in = model.add_boundary(\n",
    "    model.metabolites.amp_c,\n",
    "    boundary_type=\"demand\",\n",
    "    reaction_id=\"amp_in\")\n",
    "\n",
    "print(amp_drain)\n",
    "print(amp_in)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When a boundary reaction is created, a 'boundary metabolite' is also created as a proxy metabolite. The proxy metabolite is the external metabolite concentration (i.e., boundary condition) without instantiating a new `MassMetabolite` object to represent the external metabolite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.141257Z",
     "iopub.status.busy": "2021-01-25T18:00:06.140467Z",
     "iopub.status.idle": "2021-01-25T18:00:06.143509Z",
     "shell.execute_reply": "2021-01-25T18:00:06.143964Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'amp_b'"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "amp_in.boundary_metabolite"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The value of the 'boundary metabolite' can be set using the `MassModel.add_boundary_conditions()` method. Boundary conditions are accessed through the `MassModel.boundary_conditions` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.148886Z",
     "iopub.status.busy": "2021-01-25T18:00:06.148069Z",
     "iopub.status.idle": "2021-01-25T18:00:06.150893Z",
     "shell.execute_reply": "2021-01-25T18:00:06.151339Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'amp_b': 1.0}"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.add_boundary_conditions({amp_in.boundary_metabolite: 1})\n",
    "model.boundary_conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The automatic generation of the boundary reaction can be useful. However, sometimes the reaction stoichiometry needs to be switched in order to be intuitive. In this case, the stoichiometry of the AMP source reaction should be reversed to show that AMP enters the system, which is accomplished by using the `MassReaction.reverse_stoichiometry()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.158004Z",
     "iopub.status.busy": "2021-01-25T18:00:06.157344Z",
     "iopub.status.idle": "2021-01-25T18:00:06.159870Z",
     "shell.execute_reply": "2021-01-25T18:00:06.160318Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "amp_in:  --> amp_c\n"
     ]
    }
   ],
   "source": [
    "amp_in.reverse_stoichiometry(inplace=True)\n",
    "print(amp_in)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the addition of these two reactions adds an another pathway to the null space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.173220Z",
     "iopub.status.busy": "2021-01-25T18:00:06.171779Z",
     "iopub.status.idle": "2021-01-25T18:00:06.177145Z",
     "shell.execute_reply": "2021-01-25T18:00:06.177762Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Path 1</th>\n",
       "      <th>Path 2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>distr</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>use</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>form</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>buffer</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>amp_drain</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>amp_in</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           Path 1  Path 2\n",
       "distr           0       0\n",
       "use             1       0\n",
       "form            1       0\n",
       "buffer          0       0\n",
       "amp_drain       0       1\n",
       "amp_in          0       1"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ns = nullspace(model.S).T\n",
    "for i, row in enumerate(ns):\n",
    "    # Divide by the minimum to get all integers\n",
    "    ns[i] = np.rint(row / np.min(row[np.nonzero(row)]))\n",
    "ns = ns.T\n",
    "pd.DataFrame(\n",
    "    ns, index=model.reactions.list_attr(\"id\"),  # Rows represent reactions\n",
    "    columns=[\"Path 1\", \"Path 2\"], dtype=np.int_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining custom rates\n",
    "\n",
    "In this model, the rate for the AMP source reaction should remain at a fixed input value. However, the current rate expression for the AMP source reaction is dependent on an external AMP metabolite that exists as a boundary condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.183965Z",
     "iopub.status.busy": "2021-01-25T18:00:06.183222Z",
     "iopub.status.idle": "2021-01-25T18:00:06.185923Z",
     "shell.execute_reply": "2021-01-25T18:00:06.186516Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "amp_b*kf_amp_in\n"
     ]
    }
   ],
   "source": [
    "print(amp_in.rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, the rate can be set as a fixed input by using a custom rate expression. Custom rate expressions can be set for reactions in a model using the `MassModel.add_custom_rate()` method as follows: by passing the reaction object, a string representation of the custom rate expression, and a dictionary containing any custom parameter associated with the rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.208970Z",
     "iopub.status.busy": "2021-01-25T18:00:06.206881Z",
     "iopub.status.idle": "2021-01-25T18:00:06.211897Z",
     "shell.execute_reply": "2021-01-25T18:00:06.212424Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b1\n"
     ]
    }
   ],
   "source": [
    "model.add_custom_rate(amp_in, custom_rate=\"b1\",\n",
    "                      custom_parameters={\"b1\": 0.03})\n",
    "print(model.rates[amp_in])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ensuring Model Completeness\n",
    "### Inspecting rates and ODEs\n",
    "\n",
    "According to the [network schema](#Constructing-Models) at the start of the notebook, the network has been fully reconstructed. The reaction rates and metabolite ODEs can be inspected to ensure that the model was built without any issues.\n",
    "\n",
    "The `MassModel.rates` property is used to return a dictionary containing reactions and symbolic expressions of their rates. The model always prioritizes custom rate expressions over automatically generated mass action rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.232218Z",
     "iopub.status.busy": "2021-01-25T18:00:06.231406Z",
     "iopub.status.idle": "2021-01-25T18:00:06.234832Z",
     "shell.execute_reply": "2021-01-25T18:00:06.235386Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "distr: kf_distr*(adp_c(t)**2 - amp_c(t)*atp_c(t)/Keq_distr)\n",
      "use: kf_use*atp_c(t)\n",
      "form: kf_form*adp_c(t)\n",
      "buffer: kf_buffer*(B(t)*atp_c(t) - BP(t)*adp_c(t)/Keq_buffer)\n",
      "amp_drain: kf_amp_drain*amp_c(t)\n",
      "amp_in: b1\n"
     ]
    }
   ],
   "source": [
    "for reaction, rate in model.rates.items():\n",
    "    print(\"{0}: {1}\".format(reaction.id, rate))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, the model can access the ODEs for metabolites using the `ordinary_differential_equations` property (alias `odes`) to return a dictionary of metabolites and symbolic expressions of their ODEs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.275035Z",
     "iopub.status.busy": "2021-01-25T18:00:06.273893Z",
     "iopub.status.idle": "2021-01-25T18:00:06.278351Z",
     "shell.execute_reply": "2021-01-25T18:00:06.278799Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "adp_c: kf_buffer*(B(t)*atp_c(t) - BP(t)*adp_c(t)/Keq_buffer) - 2*kf_distr*(adp_c(t)**2 - amp_c(t)*atp_c(t)/Keq_distr) - kf_form*adp_c(t) + kf_use*atp_c(t)\n",
      "atp_c: -kf_buffer*(B(t)*atp_c(t) - BP(t)*adp_c(t)/Keq_buffer) + kf_distr*(adp_c(t)**2 - amp_c(t)*atp_c(t)/Keq_distr) + kf_form*adp_c(t) - kf_use*atp_c(t)\n",
      "amp_c: b1 - kf_amp_drain*amp_c(t) + kf_distr*(adp_c(t)**2 - amp_c(t)*atp_c(t)/Keq_distr)\n",
      "B: -kf_buffer*(B(t)*atp_c(t) - BP(t)*adp_c(t)/Keq_buffer)\n",
      "BP: kf_buffer*(B(t)*atp_c(t) - BP(t)*adp_c(t)/Keq_buffer)\n"
     ]
    }
   ],
   "source": [
    "for metabolite, ode in model.odes.items():\n",
    "    print(\"{0}: {1}\".format(metabolite.id, ode))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compartments\n",
    "\n",
    "Compartments, defined in metabolites, are recognized by the model and can be viewed in the `compartments` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.283376Z",
     "iopub.status.busy": "2021-01-25T18:00:06.282663Z",
     "iopub.status.idle": "2021-01-25T18:00:06.285506Z",
     "shell.execute_reply": "2021-01-25T18:00:06.285950Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'c': ''}"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.compartments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this model, \"c\" is an abbreviation for \"compartment\". The `compartments` attribute can be updated to reflect this mapping using a `dict`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.291894Z",
     "iopub.status.busy": "2021-01-25T18:00:06.290807Z",
     "iopub.status.idle": "2021-01-25T18:00:06.294412Z",
     "shell.execute_reply": "2021-01-25T18:00:06.294886Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'c': 'compartment'}"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.compartments = {\"c\": \"compartment\"}\n",
    "model.compartments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Units\n",
    "`Unit` and `UnitDefinition` objects are implemented as per the [SBML Unit](https://synonym.caltech.edu/software/libsbml/5.18.0/docs/formatted/python-api/classlibsbml_1_1_unit.html) and [SBML UnitDefinition](https://synonym.caltech.edu/software/libsbml/5.18.0/docs/formatted/python-api/classlibsbml_1_1_unit_definition.html) specifications.\n",
    "It can be useful for comparative reasons to create `Unit` and `UnitDefinition` objects for the model (e.g., amount, volume, time) to provide additional context. However, the model does not maintain unit consistency automatically. It is the responsibility of the users to ensure consistency among units and associated numerical values in a model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.299619Z",
     "iopub.status.busy": "2021-01-25T18:00:06.298376Z",
     "iopub.status.idle": "2021-01-25T18:00:06.300902Z",
     "shell.execute_reply": "2021-01-25T18:00:06.301371Z"
    }
   },
   "outputs": [],
   "source": [
    "from mass import Unit, UnitDefinition\n",
    "from mass.core.units import print_defined_unit_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SBML defines units using a compositional approach. The `Unit` objects represent references to base units. A `Unit` has four attributes: `kind`, `exponent`, `scale`, and `multiplier`. The `kind` attribute indicates the base unit. Valid base units are viewed using the `print_defined_unit_values()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.309117Z",
     "iopub.status.busy": "2021-01-25T18:00:06.308272Z",
     "iopub.status.idle": "2021-01-25T18:00:06.311118Z",
     "shell.execute_reply": "2021-01-25T18:00:06.311590Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "╒═════════════════════════════╕\n",
      "│ SBML Base Unit Kinds        │\n",
      "╞═════════════════════════════╡\n",
      "│ Base Unit        SBML Value │\n",
      "│ -------------  ------------ │\n",
      "│ ampere                    0 │\n",
      "│ avogadro                  1 │\n",
      "│ becquerel                 2 │\n",
      "│ candela                   3 │\n",
      "│ coulomb                   5 │\n",
      "│ dimensionless             6 │\n",
      "│ farad                     7 │\n",
      "│ gram                      8 │\n",
      "│ gray                      9 │\n",
      "│ henry                    10 │\n",
      "│ hertz                    11 │\n",
      "│ item                     12 │\n",
      "│ joule                    13 │\n",
      "│ katal                    14 │\n",
      "│ kelvin                   15 │\n",
      "│ kilogram                 16 │\n",
      "│ liter                    17 │\n",
      "│ litre                    18 │\n",
      "│ lumen                    19 │\n",
      "│ lux                      20 │\n",
      "│ meter                    21 │\n",
      "│ metre                    22 │\n",
      "│ mole                     23 │\n",
      "│ newton                   24 │\n",
      "│ ohm                      25 │\n",
      "│ pascal                   26 │\n",
      "│ radian                   27 │\n",
      "│ second                   28 │\n",
      "│ siemens                  29 │\n",
      "│ sievert                  30 │\n",
      "│ steradian                31 │\n",
      "│ tesla                    32 │\n",
      "│ volt                     33 │\n",
      "│ watt                     34 │\n",
      "│ weber                    35 │\n",
      "│ invalid                  36 │\n",
      "╘═════════════════════════════╛\n"
     ]
    }
   ],
   "source": [
    "print_defined_unit_values(\"BaseUnitKinds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `exponent`, `scale` and `multiplier` attributes indicate how the base unit should be transformed. For this model, the unit for concentration, \"Millimolar\", which is represented as millimole per liter and composed of the following base units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.317103Z",
     "iopub.status.busy": "2021-01-25T18:00:06.316198Z",
     "iopub.status.idle": "2021-01-25T18:00:06.318147Z",
     "shell.execute_reply": "2021-01-25T18:00:06.318578Z"
    }
   },
   "outputs": [],
   "source": [
    "millimole = Unit(kind=\"mole\", exponent=1, scale=-3, multiplier=1)\n",
    "per_liter = Unit(kind=\"liter\", exponent=-1, scale=1, multiplier=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combinations of `Unit` objects are contained inside a `UnitDefintion` object. `UnitDefinition` objects have three attributes: an `id`, an optional `name` to represent the combination, and a `list_of_units` attribute that contain references to the `Unit` objects. The concentration unit “Millimolar” is abbreviated as “mM” and defined as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.325612Z",
     "iopub.status.busy": "2021-01-25T18:00:06.323974Z",
     "iopub.status.idle": "2021-01-25T18:00:06.329117Z",
     "shell.execute_reply": "2021-01-25T18:00:06.329883Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Millimolar:\n",
      "<Unit at 0x7f835967d1c0 kind: mole; exponent: 1; scale: -3; multiplier: 1>\n",
      "<Unit at 0x7f835967dbe0 kind: liter; exponent: -1; scale: 1; multiplier: 1>\n"
     ]
    }
   ],
   "source": [
    "concentration_unit = UnitDefinition(id=\"mM\", name=\"Millimolar\",\n",
    "                                    list_of_units=[millimole, per_liter])\n",
    "print(\"{0}:\\n{1!r}\\n{2!r}\".format(\n",
    "    concentration_unit.name, *concentration_unit.list_of_units))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`UnitDefinition` objects also have the `UnitDefinition.create_unit()` method to directly create `Unit` objects within the `UnitDefintion`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.337986Z",
     "iopub.status.busy": "2021-01-25T18:00:06.336887Z",
     "iopub.status.idle": "2021-01-25T18:00:06.341050Z",
     "shell.execute_reply": "2021-01-25T18:00:06.341563Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min\n",
      "[<Unit at 0x7f8309bbea90 kind: second; exponent: 1; scale: 1; multiplier: 60>]\n"
     ]
    }
   ],
   "source": [
    "time_unit = UnitDefinition(id=\"min\", name=\"Minute\")\n",
    "time_unit.create_unit(kind=\"second\", exponent=1, scale=1, multiplier=60)\n",
    "print(time_unit)\n",
    "print(time_unit.list_of_units)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once created, `UnitDefintion` objects are added to the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.346625Z",
     "iopub.status.busy": "2021-01-25T18:00:06.345775Z",
     "iopub.status.idle": "2021-01-25T18:00:06.349218Z",
     "shell.execute_reply": "2021-01-25T18:00:06.348755Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<UnitDefinition Millimolar \"mM\" at 0x7f8318e6b970>,\n",
       " <UnitDefinition Minute \"min\" at 0x7f8309bbe8b0>]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.add_units([concentration_unit, time_unit])\n",
    "model.units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checking model completeness\n",
    "\n",
    "Once constructed, the model should be checked for completeness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.353906Z",
     "iopub.status.busy": "2021-01-25T18:00:06.352851Z",
     "iopub.status.idle": "2021-01-25T18:00:06.355432Z",
     "shell.execute_reply": "2021-01-25T18:00:06.355915Z"
    }
   },
   "outputs": [],
   "source": [
    "from mass import qcqa_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `qcqa_model()` function can be used to print a report about the model's completeness based on the set kwargs. The `qcqa_model()` function is used to ensure that all numerical values necessary for simulating the model are defined by setting the `parameters` and `concentrations` kwargs as `True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.380535Z",
     "iopub.status.busy": "2021-01-25T18:00:06.378942Z",
     "iopub.status.idle": "2021-01-25T18:00:06.382387Z",
     "shell.execute_reply": "2021-01-25T18:00:06.382922Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "╒══════════════════════════════════════════════╕\n",
      "│ MODEL ID: Phosphate_Trafficking              │\n",
      "│ SIMULATABLE: False                           │\n",
      "│ PARAMETERS NUMERICALY CONSISTENT: True       │\n",
      "╞══════════════════════════════════════════════╡\n",
      "│ ============================================ │\n",
      "│             MISSING PARAMETERS               │\n",
      "│ ============================================ │\n",
      "│ Reaction Parameters                          │\n",
      "│ ---------------------                        │\n",
      "│ amp_drain: kf                                │\n",
      "│ ============================================ │\n",
      "╘══════════════════════════════════════════════╛\n"
     ]
    }
   ],
   "source": [
    "qcqa_model(model, parameters=True, concentrations=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown in the report above, the forward rate constant for the AMP drain reaction was never defined. Therefore, the forward rate constant is defined, and the model is checked again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-25T18:00:06.388159Z",
     "iopub.status.busy": "2021-01-25T18:00:06.387384Z",
     "iopub.status.idle": "2021-01-25T18:00:06.390046Z",
     "shell.execute_reply": "2021-01-25T18:00:06.390492Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "╒══════════════════════════════════════════╕\n",
      "│ MODEL ID: Phosphate_Trafficking          │\n",
      "│ SIMULATABLE: True                        │\n",
      "│ PARAMETERS NUMERICALY CONSISTENT: True   │\n",
      "╞══════════════════════════════════════════╡\n",
      "╘══════════════════════════════════════════╛\n"
     ]
    }
   ],
   "source": [
    "amp_drain.kf = 0.03\n",
    "\n",
    "qcqa_model(model, parameters=True, concentrations=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the report shows that the model is not missing any values necessary for simulation. See [Checking Model Quality](./quality_assurance.ipynb) for more information on quality assurance functions and the `qcqa` submodule."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional Examples\n",
    "For additional examples on constructing models, see the following: \n",
    "\n",
    "* [Constructing Glycolysis](../gallery/workflows/constructing_glycolysis.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$^1$ Trafficking model is created from Chapter 8 of <cite data-cite=\"Pal11\">Systems Biology: Simulation of Dynamic Network States</cite>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}