{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variability Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flux Variability Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Flux variability analysis](https://pubmed.ncbi.nlm.nih.gov/20920235/) (FVA) can be used to assess the possible range of metabolic fluxes within a metabolic network while still achieving a certain level of optimality in the objective function. It essentially determines the flexibility of the metabolic network conditioned by the network's objective. The etfba package allows for FVA along with additional constraints on enzyme protein allocation and thermodynamics, implemented separately or together. The complete form (ETFVA) can be expressed as:\n",
    "<div style=\"text-align: center\">\n",
    "  <img src=\"images/ETFVA.png\" />\n",
    "</div>\n",
    "Here, $\\gamma$ in the range [0, 1] determines the fraction of the objective's optimum that should be achieved. Try a smaller $\\gamma$ if the solver cannot find a feasible solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's begin with the basic flux variability analysis, focusing solely on the mass balance constraints of metabolites. You can perform FVA using the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "model iML1515 with 2712 reactions and 1877 metabolites\n"
     ]
    }
   ],
   "source": [
    "from etfba import Model\n",
    "\n",
    "model_file = '../../models/e_coli/etfba_iML1515.bin'\n",
    "model = Model.load(model_file)\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "objective = {'BIOMASS_Ec_iML1515_core_75p37M': 1}\n",
    "flux_bound = (0, 1000)\n",
    "spec_flux_bound = {'ATPM': (6.86, 1000)}\n",
    "preset_flux = {'EX_glc__D_e_b': 10, 'FHL': 0}\n",
    "\n",
    "res = model.evaluate_variability(\n",
    "    'fva', \n",
    "    objective=objective,\n",
    "    obj_value=0.877,\n",
    "    gamma=0.99,\n",
    "    flux_bound=flux_bound,\n",
    "    spec_flux_bound=spec_flux_bound,\n",
    "    preset_flux=preset_flux\n",
    ").solve(solver='gurobi', n_jobs=36)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>Note:</b> <br></br> It is highly recommended to run variability analysis in parallel jobs by specifying the \"n_jobs\" argument, especially when dealing with large-scale models..\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated feasible range of fluxes can be accessed with the `flux_ranges` attribute, which can be further saved directly using its `save` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    CYTDK2 [ 0.000  2.331]\n",
      "      XPPT [ 0.000  1.165]\n",
      "     HXPRT [ 0.000  1.165]\n",
      "     NDPK5 [-1.576  1.186]\n",
      "    SHK3Dr [ 0.331  0.388]\n",
      "     NDPK6 [ 0.000  0.643]\n",
      "     NDPK8 [-1.576  1.185]\n",
      "    DHORTS [-0.445 -0.287]\n",
      "     OMPDC [ 0.287  0.445]\n",
      "    PYNP2r [-2.331  1.188]\n"
     ]
    }
   ],
   "source": [
    "for rxnid, flux_range in list(res.flux_ranges.items())[:10]:\n",
    "    print(f'{rxnid:>10} [{flux_range[0]:>6.3f} {flux_range[1]:>6.3f}]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we can extend the analysis to include constraints on enzyme protein allocation and thermodynamics, namely EFVA and TFVA, respectively.\n",
    "Here's how you can perform EFVA:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "kcat_file = '../../models/e_coli/kcats.xlsx'\n",
    "mw_file = '../../models/e_coli/mws.xlsx'\n",
    "dgpm_file = '../../models/e_coli/dgpms.xlsx'\n",
    "kcats = pd.read_excel(kcat_file, header=None, index_col=0).squeeze()\n",
    "mws = pd.read_excel(mw_file, header=None, index_col=0).squeeze()\n",
    "dgpms = pd.read_excel(dgpm_file, header=None, index_col=0).squeeze()\n",
    "\n",
    "# set reactions with available kcats and MWs\n",
    "eff_rxns = []\n",
    "for rxnid, rxn in model.reactions.items():\n",
    "    if rxn.rev:\n",
    "        if rxnid+'_f' in kcats.index and rxnid+'_b' in kcats.index and rxnid in mws.index:\n",
    "            rxn.forward_kcat = kcats[rxnid+'_f']\n",
    "            rxn.backward_kcat = kcats[rxnid+'_b']\n",
    "            rxn.molecular_weight = mws[rxnid]\n",
    "            eff_rxns.append(rxnid)\n",
    "    else:\n",
    "        if rxnid in kcats.index and rxnid in mws.index:\n",
    "            rxn.forward_kcat = kcats[rxnid]\n",
    "            eff_rxns.append(rxnid)\n",
    "enz_ub = 0.15\n",
    "\n",
    "# set reactions with available ΔG'm\n",
    "ex_rxns = []\n",
    "for rxnid, rxn in model.reactions.items():\n",
    "    if rxnid in dgpms.index:\n",
    "        rxn.standard_gibbs_energy = dgpms[rxnid]\n",
    "    else:\n",
    "        ex_rxns.append(rxnid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform EFVA\n",
    "preset_flux = {'FHL': 0}\n",
    "res = model.evaluate_variability(\n",
    "    'efva', \n",
    "    objective=objective,\n",
    "    obj_value=0.866,\n",
    "    gamma=0.99,\n",
    "    flux_bound=flux_bound, \n",
    "    spec_flux_bound=spec_flux_bound,\n",
    "    preset_flux=preset_flux,\n",
    "    inc_enz_cons=eff_rxns,\n",
    "    enz_prot_lb=enz_ub\n",
    ").solve(solver='gurobi', n_jobs=36)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And TFVA:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform TFVA\n",
    "preset_flux = {'EX_glc__D_e_b': 10, 'FHL': 0}\n",
    "conc_bound = (0.0001, 100)\n",
    "spec_conc_bound = {\n",
    "    'o2_c': (0.0001, 0.0082),   \n",
    "    'co2_c': (0.1, 100)   \n",
    "}\n",
    "preset_conc = {\n",
    "    'glc__D_p': 20,\n",
    "    'pi_p': 56,\n",
    "    'so4_p': 3,\n",
    "    'nh4_p': 19,\n",
    "    'na1_p': 160,\n",
    "    'k_p': 22,\n",
    "    'fe2_p': 62\n",
    "}\n",
    "res = model.evaluate_variability(\n",
    "    'tfva', \n",
    "    objective=objective,\n",
    "    obj_value=0.877,\n",
    "    gamma=0.99,\n",
    "    flux_bound=flux_bound,\n",
    "    conc_bound=conc_bound,\n",
    "    spec_flux_bound=spec_flux_bound,\n",
    "    spec_conc_bound=spec_conc_bound,\n",
    "    preset_flux=preset_flux,\n",
    "    preset_conc=preset_conc,\n",
    "    ex_thermo_cons=ex_rxns\n",
    ").solve(solver='gurobi', n_jobs=36)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moreover, we can conduct flux variability analysis with comprehensive constraints on enzyme protein allocation and thermodynamics, i.e., ETFVA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform ETFVA\n",
    "preset_flux = {'FHL': 0}\n",
    "res = model.evaluate_variability(\n",
    "    'etfva', \n",
    "    objective=objective,\n",
    "    obj_value=0.864,\n",
    "    gamma=0.99,\n",
    "    flux_bound=flux_bound,\n",
    "    conc_bound=conc_bound,\n",
    "    spec_flux_bound=spec_flux_bound,\n",
    "    spec_conc_bound=spec_conc_bound,\n",
    "    preset_flux=preset_flux,\n",
    "    preset_conc=preset_conc,\n",
    "    ex_thermo_cons=ex_rxns,\n",
    "    inc_enz_cons=eff_rxns,\n",
    "    enz_prot_lb=enz_ub\n",
    ").solve(solver='gurobi', n_jobs=36)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enzyme Protein Cost Variability Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to flux variability analysis, the variability of enzyme protein costs in enzymatic reactions can also be evaluated. Below is the complete form of enzyme protein variability analysis (TEVA), where the constraints are identical to those in ETFVA.\n",
    "<div style=\"text-align: center\">\n",
    "  <img src=\"images/TEVA.gif\" />\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If only the constraint of enzyme protein is considered, basic EVA can be performed using the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "preset_flux = {'FHL': 0}\n",
    "res = model.evaluate_variability(\n",
    "    'eva', \n",
    "    objective=objective,\n",
    "    obj_value=0.866,\n",
    "    gamma=0.99,\n",
    "    flux_bound=flux_bound, \n",
    "    spec_flux_bound=spec_flux_bound,\n",
    "    preset_flux=preset_flux,\n",
    "    inc_enz_cons=eff_rxns,\n",
    "    enz_prot_lb=enz_ub\n",
    ").solve(solver='gurobi', n_jobs=36)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated feasible ranges of enzyme protein costs can be found in the `protein_cost_ranges` attribute. You can use the `save` method of this attribute to save the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       ENO [0.00875 0.01220]\n",
      "      GAPD [0.00612 0.00885]\n",
      "     KARA2 [0.00499 0.00641]\n",
      "     KARA1 [0.00490 0.00632]\n",
      "       PGK [0.00489 0.00736]\n",
      "      METS [0.00474 0.00548]\n",
      "       FBA [0.00357 0.01068]\n",
      "       PGM [0.00289 0.00498]\n",
      "    GHMT2r [0.00249 0.00504]\n",
      "  ATPS4rpp [0.00196 0.00346]\n"
     ]
    }
   ],
   "source": [
    "pro_cost_sorted = sorted(\n",
    "    res.protein_cost_ranges.items(), \n",
    "    key=lambda item: item[1][0], \n",
    "    reverse=True\n",
    ")\n",
    "for rxnid, pro_cost_range in pro_cost_sorted[:10]:\n",
    "    print(f'{rxnid:>10} [{pro_cost_range[0]:>6.5f} {pro_cost_range[1]:>6.5f}]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below shows the code of enzyme protein variability analysis under the complete constraints of enzyme protein allocation and thermodynamics, i.e., TEVA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = model.evaluate_variability(\n",
    "    'teva', \n",
    "    objective=objective,\n",
    "    obj_value=0.864,\n",
    "    gamma=0.98,\n",
    "    flux_bound=flux_bound,\n",
    "    conc_bound=conc_bound,\n",
    "    spec_flux_bound=spec_flux_bound,\n",
    "    spec_conc_bound=spec_conc_bound,\n",
    "    preset_flux=preset_flux,\n",
    "    preset_conc=preset_conc,\n",
    "    ex_thermo_cons=ex_rxns,\n",
    "    inc_enz_cons=eff_rxns,\n",
    "    enz_prot_lb=enz_ub\n",
    ").solve(solver='gurobi', n_jobs=36)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thermodynamic Variability Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thermodynamic variability analysis can be used to evaluate the feasible ranges of reaction Gibbs energy change. Its complete form (ETVA) can be represented as below, with constraints identical to those in ETFVA.\n",
    "<div style=\"text-align: center\">\n",
    "  <img src=\"images/ETVA.gif\" />\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basic version of TVA without enzyme protein allocation constraint can be conducted with the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "preset_flux = {'EX_glc__D_e_b': 10, 'FHL': 0}\n",
    "res = model.evaluate_variability(\n",
    "    'tva',\n",
    "    objective=objective,\n",
    "    obj_value=0.877,\n",
    "    gamma=0.99,\n",
    "    flux_bound=flux_bound,\n",
    "    conc_bound=conc_bound,\n",
    "    spec_flux_bound=spec_flux_bound,\n",
    "    spec_conc_bound=spec_conc_bound,\n",
    "    preset_flux=preset_flux,\n",
    "    preset_conc=preset_conc,\n",
    "    ex_thermo_cons=ex_rxns\n",
    ").solve(solver='gurobi', n_jobs=36)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To inspect the bounds of feasible $\\Delta G'$, one can use the `gibbs_energy_ranges` attribute, which also has a `save` method for direct saving."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    CYTDK2 [-102.0   44.3]\n",
      "      XPPT [ -93.9   34.7]\n",
      "     HXPRT [ -97.7   31.0]\n",
      "     NDPK5 [ -60.0   41.8]\n",
      "    SHK3Dr [ -56.7   -0.0]\n",
      "     NDPK6 [ -63.4   38.4]\n",
      "     NDPK8 [ -59.6   42.3]\n",
      "    DHORTS [   0.0   17.1]\n",
      "     OMPDC [ -69.4   -0.0]\n",
      "    PYNP2r [ -68.5   67.5]\n"
     ]
    }
   ],
   "source": [
    "for rxnid, dgp_range in list(res.gibbs_energy_ranges.items())[:10]:\n",
    "    print(f'{rxnid:>10} [{dgp_range[0]:>6.1f} {dgp_range[1]:>6.1f}]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Variability of reactions Gibbs energy under complete constraints of enzyme protein allocation and thermodynamics, i.e., ETVA, can be conducted as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "preset_flux = {'FHL': 0}\n",
    "res = model.evaluate_variability(\n",
    "    'etva', \n",
    "    objective=objective,\n",
    "    obj_value=0.864,\n",
    "    gamma=0.98,\n",
    "    flux_bound=flux_bound,\n",
    "    conc_bound=conc_bound,\n",
    "    spec_flux_bound=spec_flux_bound,\n",
    "    spec_conc_bound=spec_conc_bound,\n",
    "    preset_flux=preset_flux,\n",
    "    preset_conc=preset_conc,\n",
    "    ex_thermo_cons=ex_rxns,\n",
    "    inc_enz_cons=eff_rxns,\n",
    "    enz_prot_lb=enz_ub\n",
    ").solve(solver='gurobi', n_jobs=36)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "etfba-py38",
   "language": "python",
   "name": "etfba-py38"
  },
  "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.8.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}