{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building a Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building with Metabolite and Reaction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Building a model starts with creating metabolites and reactions. Let’s consider the reaction OAA(abcd) + AcCoA(ef) $\\rightarrow$ Cit(dcbfea) in the toy model as an example. To create the reactants, we can use the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Metabolite OAA(a,b,c,d) Metabolite AcCoA(d,f) Metabolite Cit(d,c,b,f,e,a)\n"
     ]
    }
   ],
   "source": [
    "from freeflux import Metabolite, Reaction, Model\n",
    "\n",
    "oaa = Metabolite('OAA', atoms = ['a', 'b', 'c', 'd'])\n",
    "accoa = Metabolite('AcCoA', atoms = ['d', 'f'])\n",
    "cit = Metabolite('Cit', atoms = list('dcbfea'))\n",
    "print(oaa, accoa, cit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "<b>Note:</b> <br></br> Rotationally symmetric metabolites usually have <a href=\"https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1994654/\">equivalents</a> in biochemical reactions, e.g., succinate and fumarate in the TCA cycle. In this case, the \"atoms\" argument in Metabolite constructor needs to be assigned as \"abcd,dcba\". \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we can build a reaction that consumes and produces these metabolites:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reaction v1: AcCoA+OAA->Cit\n"
     ]
    }
   ],
   "source": [
    "v1 = Reaction('v1', reversible = False)\n",
    "v1.add_substrates([oaa, accoa], stoichiometry = [1, 1])\n",
    "v1.add_products(cit, stoichiometry = 1)\n",
    "print(v1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we add the reaction to the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model demo (3 metabolites, 1)\n",
      "Reaction v1: AcCoA+OAA->Cit\n"
     ]
    }
   ],
   "source": [
    "demo = Model('demo')\n",
    "demo.add_reactions([v1])\n",
    "print(demo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The process can be repeated until all reactions are involved."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To modify the mode, one can use `remove_substrates` and `remove_products` to delete reactants from a reaction, and `remove_reactions` to delate reactions from a model. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading from File"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To input a large metabolic network, it is convenient to load from a file. The following formats are supported: tab-separated values (.tsv) or Excel spreadsheet (.xlsx). The file should have the following format:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|#reaction_ID|reactant_IDs(atom)|product_IDs(atom)|reversibility|\n",
    "|:------:|:------:|:------:|:------:|\n",
    "|v1|OAA(abcd)+AcCoA(ef)|Cit(dcbfea)|0|\n",
    "|v2|Cit(abcdef)|AKG(abcde)+CO2(f)|0|\n",
    "|v3|AKG(abcde)|Glu(abcde)|0|\n",
    "|v4|AKG(abcde)|Suc(bcde)+CO2(a)|0|\n",
    "|v5|Suc(abcd,dcba)|Fum(abcd,dcba)|0|\n",
    "|v6|Fum(abcd,dcba)|OAA(abcd)|1|\n",
    "|v7|Asp(abcd)|OAA(abcd)|0|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "<b>Note:</b> <br></br> 1. \"#\" is required in the header; <br></br> 2. Metabolite name could include but must not start with digits; <br></br> 3. Reactions with end metabolite, i.e., initial substrates and final products of the network (for example, substrate glucose and biomass) should be irreversible.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To load the model from the file, use the `read_from_file` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model demo (9 metabolites, 7)\n",
      "Reaction v1: AcCoA+OAA->Cit\n",
      "Reaction v2: Cit->AKG+CO2\n",
      "Reaction v3: AKG->Glu\n",
      "Reaction v4: AKG->CO2+Suc\n",
      "Reaction v5: Suc->Fum\n",
      "Reaction v6: Fum<->OAA\n",
      "Reaction v7: Asp->OAA\n"
     ]
    }
   ],
   "source": [
    "MODEL_FILE = 'path/to/reactions.tsv'\n",
    "\n",
    "demo.read_from_file(MODEL_FILE)\n",
    "print(demo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model file can be found [here](https://github.com/Chaowu88/freeflux/tree/main/models/toy)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `metabolites_info` and `reactions_info` attributes can be used to access information about the metabolites and reactions in the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'OAA': [Metabolite OAA(abcd)],\n",
       " 'AcCoA': [Metabolite AcCoA(ef)],\n",
       " 'Cit': [Metabolite Cit(dcbfea), Metabolite Cit(abcdef)],\n",
       " 'AKG': [Metabolite AKG(abcde)],\n",
       " 'CO2': [Metabolite CO2(f), Metabolite CO2(a)],\n",
       " 'Glu': [Metabolite Glu(abcde)],\n",
       " 'Suc': [Metabolite Suc(bcde), Metabolite Suc(abcd,dcba)],\n",
       " 'Fum': [Metabolite Fum(abcd,dcba)],\n",
       " 'Asp': [Metabolite Asp(abcd)]}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "demo.metabolites_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OrderedDict([('v1', Reaction v1: AcCoA+OAA->Cit),\n",
       "             ('v2', Reaction v2: Cit->AKG+CO2),\n",
       "             ('v3', Reaction v3: AKG->Glu),\n",
       "             ('v4', Reaction v4: AKG->CO2+Suc),\n",
       "             ('v5', Reaction v5: Suc->Fum),\n",
       "             ('v6', Reaction v6: Fum<->OAA),\n",
       "             ('v7', Reaction v7: Asp->OAA)])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "demo.reactions_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain the stoichiometric matrix for the net reactions, use the `get_net_stoichiometric_matrix` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "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>v1</th>\n",
       "      <th>v2</th>\n",
       "      <th>v3</th>\n",
       "      <th>v4</th>\n",
       "      <th>v5</th>\n",
       "      <th>v6</th>\n",
       "      <th>v7</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AKG</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Cit</th>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Fum</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>OAA</th>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Suc</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      v1   v2   v3   v4   v5   v6   v7\n",
       "AKG  0.0  1.0 -1.0 -1.0  0.0  0.0  0.0\n",
       "Cit  1.0 -1.0  0.0  0.0  0.0  0.0  0.0\n",
       "Fum  0.0  0.0  0.0  0.0  1.0 -1.0  0.0\n",
       "OAA -1.0  0.0  0.0  0.0  0.0  1.0  1.0\n",
       "Suc  0.0  0.0  0.0  1.0 -1.0  0.0  0.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "demo.get_net_stoichiometric_matrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, the stoichiometric matrix for the total reactions can be obtained using the `get_total_stoichiometric_matrix` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "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>v1</th>\n",
       "      <th>v2</th>\n",
       "      <th>v3</th>\n",
       "      <th>v4</th>\n",
       "      <th>v5</th>\n",
       "      <th>v6_f</th>\n",
       "      <th>v6_b</th>\n",
       "      <th>v7</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AKG</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Cit</th>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Fum</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>OAA</th>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Suc</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      v1   v2   v3   v4   v5  v6_f  v6_b   v7\n",
       "AKG  0.0  1.0 -1.0 -1.0  0.0   0.0  -0.0  0.0\n",
       "Cit  1.0 -1.0  0.0  0.0  0.0   0.0  -0.0  0.0\n",
       "Fum  0.0  0.0  0.0  0.0  1.0  -1.0   1.0  0.0\n",
       "OAA -1.0  0.0  0.0  0.0  0.0   1.0  -1.0  1.0\n",
       "Suc  0.0  0.0  0.0  1.0 -1.0   0.0  -0.0  0.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "demo.get_total_stoichiometric_matrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metabolic Network Decomposition\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to establish the mapping between fluxes and labeling patterns of metabolites, or more specifically, the elementary metabolite units [(EMU)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1994654/), a metabolic network carrying atom transfer information needs to be decomposed. This can be achieved using the [adjacency matrix enabled method](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6503117/) that has been previously proposed and is implemented by the `decompose_network` function.\n",
    "To decompose the network, one can use the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "emu_mats = demo.decompose_network(\n",
    "    ini_emus = {'Glu': ['12345']}   # define initial EMU Glu_12345 for decomposition\n",
    ")   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function returns a dictionary of the decomposed EMU networks with the size as the key."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "<b>Note:</b> <br></br> The `decompose_network` can run with parallel jobs for multiple target EMUs by specifing the \"n_jobs\" argument. \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, to obtain the EMU adjacency matrix for EMUs of size 1, the following code can be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "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>EMU Fum_2</th>\n",
       "      <th>EMU OAA_2</th>\n",
       "      <th>EMU OAA_3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>EMU Fum_2</th>\n",
       "      <td>0</td>\n",
       "      <td>1.0*v6_f</td>\n",
       "      <td>1.0*v6_f</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>EMU OAA_2</th>\n",
       "      <td>0.5*v5 + 0.5*v6_b</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>EMU OAA_3</th>\n",
       "      <td>0.5*v6_b</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>EMU AcCoA_2</th>\n",
       "      <td>0.5*v5</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>EMU Asp_2</th>\n",
       "      <td>0</td>\n",
       "      <td>1.0*v7</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>EMU Asp_3</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0*v7</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     EMU Fum_2 EMU OAA_2 EMU OAA_3\n",
       "EMU Fum_2                    0  1.0*v6_f  1.0*v6_f\n",
       "EMU OAA_2    0.5*v5 + 0.5*v6_b         0         0\n",
       "EMU OAA_3             0.5*v6_b         0         0\n",
       "EMU AcCoA_2             0.5*v5         0         0\n",
       "EMU Asp_2                    0    1.0*v7         0\n",
       "EMU Asp_3                    0         0    1.0*v7"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "emu_mats[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This EMU adjacency matrix shows the transformation of size-1 EMUs bridged by reactions. For instance, the EMU OAA_2 comes from precursor EMU Fum_2 through the forward reaction v6_f and the substrate EMU Asp_2 through the reaction v7. Through the decomposition, the labeling pattern of any target EMU(s) can be traced back to the initial labeled or unlabeled substrate(s)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:freeflux-py37]",
   "language": "python",
   "name": "conda-env-freeflux-py37-py"
  },
  "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}