{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Getting Started\n",
    "\n",
    "[multiTFA](https://github.com/biosustain/multitfa) is a python package for accurate thermodynamic flux analysis. multiTFA is built on COBRApy to take advantage of COBRA model structure. Installation instructions can be found in **README.rst**. \n",
    "\n",
    "## 1.1. Building a multiTFA model\n",
    "\n",
    "To create a multitfa model, we require \n",
    "\n",
    "1. A COBRA model\n",
    "\n",
    "2. Metabolite database identifiers, such as BIGG, KEGG, SEED etc. - This is required to match and retrieve thermodyanamic information against our database. Metabolite identifiers can be added to the multiTFA model later. For simplicity, we use the same nomneclature as [equilibrator-api](https://gitlab.com/equilibrator/equilibrator-api). \n",
    "\n",
    "3. Compartment specific pH and ionic strength information - To calculate Gibbs energy related to transmembrane transport, information about membrane potential is also required in the form of pandas Dataframes. \n",
    "\n",
    "Optionally, intracellular metabolite concentrations can be provided to calculate the transformed Gibbs energy of reactions. If metabolite concentrations are not available, a broad metabolite concentration range is assumed.   \n",
    "\n",
    "**NOTE**: All the input information should be in S.I units, e.g: ionic strength, metabolite concentration in Moles.\n",
    "\n",
    "To demonstrate the usage, we use *E. Coli* core model from [BIGG database](http://bigg.ucsd.edu/).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First step in creating a multitfa model is to get a COBRA model. For this example, we load *E. coli* core model. This is available in `Examples` directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra.io import load_model\n",
    "\n",
    "# Load the Cobra model\n",
    "cobra_model = load_model('e_coli_core')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us define the compartment properties such as pH, ionic strength (M) etc. For transport calculations, we will provide membrane potential in mV. For convenience, we use `pandas DataFrame` to define these properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compartment specific pH and Ionic strength information. Ionic strength units should be in M.\n",
    "pH_I_dict = {\n",
    "    \"pH\": {\"c\": 7.5, \"e\": 7, \"p\": 7},\n",
    "    \"I\": {\"c\": 0.25, \"e\": 0, \"p\": 0},\n",
    "}\n",
    "\n",
    "# Membrane potentials between different compartments. Units for membrane potential is mV.\n",
    "del_psi_dict = { \n",
    "\"c\": {\"c\":0,\"e\": 0, \"p\": 150}, \n",
    "\"e\": {\"c\": 0,\"e\":0, \"p\": 0}, \n",
    "\"p\": {\"c\": -150, \"e\": 0, \"p\":0}, \n",
    "} \n",
    "\n",
    "import pandas as pd\n",
    "# Define membrane potential and compartment information in Dataframes\n",
    "del_psi = pd.DataFrame.from_dict(data=del_psi_dict)\n",
    "comp_info = pd.DataFrame.from_dict(data=pH_I_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, it is time to list any reactions that needs to be excluded from thermodynamic analysis. For example, Biomass, demand reactions etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# List of reaction ids that you would like to exclude from thermodynamic analysis, here we exclude demand reactions\n",
    "Excl = [\n",
    "    rxn.id\n",
    "    for rxn in cobra_model.reactions\n",
    "    if rxn.id.startswith(\"EX_\") or rxn.id.startswith(\"DM_\")\n",
    "] + [\"BIOMASS_Ecoli_core_w_GAM\", \"O2t\", \"H2Ot\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we have all the information we need. We can now create an instance of multitfa model.\n",
    "\n",
    "Please note, when you use multitfa for the very first time, it downloads datafiles of about 3 GB. This is a one time thing only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading package metadata...\n",
      "Fragments already downloaded\n",
      "Downloading package metadata...\n",
      "Fragments already downloaded\n",
      "Downloading package metadata...\n",
      "Fragments already downloaded\n"
     ]
    }
   ],
   "source": [
    "# Now build the multiTFA model from COBRA model\n",
    "from multitfa.core import tmodel\n",
    "tfa_model = tmodel(\n",
    "    cobra_model, Exclude_list=Excl, compartment_info=comp_info, membrane_potential=del_psi\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to access the thermodynamic information from equilibrator api, we need to add a external database identifier against all the metabolites in tfa model. This can be done by populating the 'Kegg_id' property of the metabolites. If you do not know the identifier for some metabolites, it is okay, you can simply add `NA` against it. For simplicity, we use the notation similar to equilibrator api. For the example below in the *E. coli* core model, we use bigg identifiers. \n",
    "\n",
    "You can follow equilibrator api notation for other database identifiers, for example, if you are using `KEGG` identifier, simply use `kegg:kegg identifier`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "for met in tfa_model.metabolites: \n",
    "    kegg_id = 'bigg.metabolite:'+met.id[:-2] \n",
    "    met.Kegg_id = kegg_id\n",
    "\n",
    "tfa_model.update()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have setup the model and populated the thermodynamic constraints, its better to check if constraints are added properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of constraints before update 72\n",
      "\n",
      "Number of reactions in the model 95\n",
      "\n",
      "Number of reactions excluded from thermodynamic analysis 23\n",
      "\n",
      "Number of mass balance constraints 72\n",
      "\n",
      "Number of constraints to be present after update 504\n",
      "\n",
      "Number of constraints present after update 504\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"Number of constraints before update {}\\n\".format(len(cobra_model.constraints)))\n",
    "print(\"Number of reactions in the model {}\\n\".format(len(tfa_model.reactions)))\n",
    "print(\"Number of reactions excluded from thermodynamic analysis {}\\n\".format(len(tfa_model.Exclude_reactions)))\n",
    "print(\"Number of mass balance constraints {}\\n\".format(len(tfa_model.metabolites)))\n",
    "num_cons = 2*3*(len(tfa_model.reactions) - len(tfa_model.Exclude_reactions)) + len(tfa_model.metabolites)\n",
    "print(\"Number of constraints to be present after update {}\\n\".format(num_cons))\n",
    "print(\"Number of constraints present after update {}\\n\".format(len(tfa_model.constraints)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 Metabolites\n",
    "\n",
    "Reactions and metabolites in `multiTFA` model can be accessed similar to a COBRA model. For demonstration, we will consider cytosolic `ATP` as our metabolite with `atp_c` identifier in t_model. We will access `ATP` and print it's standard Gibbs energy of formation and Also, change its intracellular concentration lower and upper bounds to '2e-3' and '5e-2' M respectively from it's default bounds. These changes will reflect in solver variables and constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Standard Gibbs energy of formation of ATP is -2259.1882733696866 KJ/mol\n",
      "initial concentrations\n",
      "atp_c\t1e-05\t0.02\n",
      "\n",
      "Selected objective sense:  MINIMIZE\n",
      "Selected objective  name:  f3d52b45-5dba-11eb-bac2-c1c11d7a5b17\n",
      "Selected RHS        name:  rhs\n",
      "Selected bound      name:  bnd\n",
      "Concentrations after modifying\n",
      "atp_c\t0.002\t0.05\n"
     ]
    }
   ],
   "source": [
    "atp = tfa_model.metabolites.get_by_id('atp_c')\n",
    "print(\"Standard Gibbs energy of formation of ATP is {} KJ/mol\".format(atp.delG_f))\n",
    "print(\"initial concentrations\")\n",
    "print('{}\\t{}\\t{}'.format(atp.id, atp.concentration_min, atp.concentration_max))\n",
    "atp.concentration_min = 2e-3\n",
    "atp.concentration_max = 5e-2\n",
    "print(\"Concentrations after modifying\")\n",
    "print('{}\\t{}\\t{}'.format(atp.id, atp.concentration_min, atp.concentration_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can always change the Gibbs energy of formation value of a metabolite by resetting the `delG_f` property of the metabolite. Also, you can access additional information like, major pseudoisomer at given compartment conditions (`major_ms`), equilibrator accession of the metabolite."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 Reactions\n",
    "\n",
    "Reactions in multiTFA model are of type cobra.DictList. Similar to COBRA model, reactions can be accessed/retrieved by its index or by it's name and properties/attributes can be accessed. Let us consider ATP synthase raction in our model by 'ATPS4rpp' and access it's attributes.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "adp_c + 4.0 h_e + pi_c <=> atp_c + h2o_c + 3.0 h_c\n"
     ]
    }
   ],
   "source": [
    "atps = tfa_model.reactions.get_by_id('ATPS4r')\n",
    "print(atps.reaction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Access and change bounds of `atps`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lower bound is -1000.0 and upper bound is 1000.0\n",
      "\n",
      "updated lower and upper bounds are 0 \t 1000.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"lower bound is {} and upper bound is {}\\n\".format(atps.lower_bound, atps.upper_bound))\n",
    "atps.lower_bound = 0\n",
    "print(\"updated lower and upper bounds are {} \\t {}\\n\".format(atps.lower_bound, atps.upper_bound))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access the Gibbs energy of transport and transformed Gibbs energy of the reaction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gibbs energy of transport is -4.605170185988092 KJ/mol\n",
      "Transform Gibbs energy of is 32.39600409402766 KJ/mol\n"
     ]
    }
   ],
   "source": [
    "print(\"Gibbs energy of transport is {} KJ/mol\".format(atps.delG_transport))\n",
    "print(\"Transform Gibbs energy of is {} KJ/mol\".format(atps.delG_prime))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can always reset the transport and transformed Gibbs energy of reactions. You can also access the stoichiometric matrix of the reaction, all the thermodynamic constraints associated with the reaction. Please refer to api for all the attributes of the model."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}