community_model_25.ipynb

{
 "metadata": {
  "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.6"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "python376jvsc74a57bd012dc29294301b76d7961fc9a20d6409e84490a3af29fdd65fcbf33ee51692f0a",
   "display_name": "Python 3.7.6 64-bit ('base': conda)"
  },
  "metadata": {
   "interpreter": {
    "hash": "a8f61be024eba58adef938c9aa1e29e02cb3dece83a5348b1a2dafd16a070453"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import cobra\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.optimize import linprog\n",
    "import scipy\n",
    "import mip\n",
    "from copy import deepcopy\n",
    "\n",
    "def save_dict(data, name):\n",
    "    with open(name, 'w' ) as file:\n",
    "        json.dump( data, file )"
   ]
  },
  {
   "source": [
    "# Community modeling\n",
    "\n",
    "In this notebook we will implement a method to create community models of two or more species specific metabolic models using cobrapy."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Growth:  0.282365392532252\nGrowth:  2.558694612613397\n"
     ]
    }
   ],
   "source": [
    "model_DP = cobra.io.read_sbml_model(\"models/consistent_DP_SNM.xml\")\n",
    "model_SA = cobra.io.read_sbml_model(\"models/consistent_iYS854_SNM.xml\")\n",
    "print(\"Growth: \", model_DP.slim_optimize())\n",
    "print(\"Growth: \", model_SA.slim_optimize())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "for rec in model_SA.reactions:\n",
    "    rec.lower_bound = max(rec.lower_bound, -1000)\n",
    "    rec.upper_bound = min(rec.upper_bound, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "        Compound    BiGG ModelSeed    KEGG\n",
       "0        Alanine  ala__L  cpd00035  C00041\n",
       "1       Arginine  arg__L  cpd00051  C00062\n",
       "2       Cysteine  cys__L  cpd00084  C00097\n",
       "3  Glutamic acid  glu__L  cpd00023  C00025\n",
       "4        Glycine     gly  cpd00033  C00037"
      ],
      "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>Compound</th>\n      <th>BiGG</th>\n      <th>ModelSeed</th>\n      <th>KEGG</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>Alanine</td>\n      <td>ala__L</td>\n      <td>cpd00035</td>\n      <td>C00041</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>Arginine</td>\n      <td>arg__L</td>\n      <td>cpd00051</td>\n      <td>C00062</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>Cysteine</td>\n      <td>cys__L</td>\n      <td>cpd00084</td>\n      <td>C00097</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>Glutamic acid</td>\n      <td>glu__L</td>\n      <td>cpd00023</td>\n      <td>C00025</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>Glycine</td>\n      <td>gly</td>\n      <td>cpd00033</td>\n      <td>C00037</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "source": [
    "snm3 = pd.read_csv(\"SNM3.csv\", sep =\";\")\n",
    "snm3.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "BIOMASS_DP = \"Growth\" \n",
    "BIOMASS_SA = \"BIOMASS_iYS_wild_type\"\n",
    "models = [model_DP.copy(), model_SA.copy()]"
   ]
  },
  {
   "source": [
    "The following metabolites have been previously added to obtain growth on the medium: EX_ile__L_e, EX_met__L_e and EX_26dap__M_e for DP. In the community this may is not necessary as the metabolites could be produced from e.g. SA.\n",
    "\n",
    "**PROBLEM:** SA does not have EX_26dap__M_e and hence cannot produce it for DP....\n",
    "\n",
    "Thus we cannot drop EX_26dap__M_e from the medium otherwise DP could not grow..."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "EX_met__L_e: met__L_e --> \nEX_ile__L_e: ile__L_e --> \n"
     ]
    }
   ],
   "source": [
    "print(model_SA.exchanges.EX_met__L_e)\n",
    "print(model_SA.exchanges.EX_ile__L_e)\n",
    "#print(model_SA.exchanges.EX_26dap__M_e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "10.0\n"
     ]
    }
   ],
   "source": [
    "# Can produces met\n",
    "with model_SA as model: \n",
    "    model.objective = model.exchanges.EX_met__L_e\n",
    "    print(model.slim_optimize())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "21.406249999999982\n"
     ]
    }
   ],
   "source": [
    "# Can produces ile\n",
    "with model_SA as model: \n",
    "    model.objective = model.exchanges.EX_ile__L_e\n",
    "    print(model.slim_optimize())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra.medium import minimal_medium\n",
    "\n",
    "min_medium_DP = minimal_medium(model_DP, 0.1*model_DP.slim_optimize())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_medium_SA = minimal_medium(model_SA, 0.1*model_SA.slim_optimize())"
   ]
  },
  {
   "source": [
    "## 1) Constructing of community model explicitely"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "Here we introduce only shuttle reactions for reactions that are common in the uptake/sekretion reactions of the individual models!"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_stoichiometry_matrix(model):\n",
    "    metabolites = model.metabolites \n",
    "    reactions = model.reactions \n",
    "    S = np.zeros((len(metabolites), len(reactions)))\n",
    "    \n",
    "    met_id = dict()\n",
    "    rec_id = dict()\n",
    "    for i,reaction in enumerate(model.reactions):\n",
    "        rec_id[reaction.id] = i\n",
    "        for metabolite, stoich in reaction.metabolites.items():\n",
    "            met_id[metabolite.id] = int(metabolites.index(metabolite))\n",
    "            S[metabolites.index(metabolite), i] = stoich\n",
    "    return S, met_id, rec_id "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Model():\n",
    "    def __init__(self, model, biomass_function):\n",
    "        \"\"\" This is a new class of metabolic model, capable of flux balance analysis\n",
    "        Attributes:\n",
    "        models (list): CobraPy models of single organisms which will be used in construction\n",
    "        biomass_reactions (list): List of strings containing the ids for the growth reactions\n",
    "        \"\"\"\n",
    "        self.biomass_function = biomass_function\n",
    "        self.model = model\n",
    "        self.id = model.id\n",
    "        # Compute stoichimetry_matrix\n",
    "        S, met_id, rec_id = create_stoichiometry_matrix(model)\n",
    "        self.num_reactions = S.shape[1]\n",
    "        self.num_metabolites = S.shape[0]\n",
    "        self.stoichiometry_matrix = scipy.sparse.csr_matrix(S)\n",
    "        self.met_id = met_id\n",
    "        self.rec_id = rec_id \n",
    "        # Set objective\n",
    "        idx = self.rec_id[biomass_function]\n",
    "        c = np.zeros(self.num_reactions)\n",
    "        c[idx] = 1\n",
    "        self.objective_c = c\n",
    "        # Set bounds\n",
    "        self._reset_bounds()\n",
    "    \n",
    "    @property\n",
    "    def reactions(self):\n",
    "        return self.model.reactions\n",
    "    @property\n",
    "    def exchanges(self):\n",
    "        return self.model.exchanges\n",
    "    @property\n",
    "    def metabolites(self):\n",
    "        return self.model.metabolites\n",
    "    @property\n",
    "    def medium(self):\n",
    "        return self.model.medium\n",
    "\n",
    "    def set_medium(self, medium):\n",
    "        ex_ids = [ex.id for ex in self.exchanges]\n",
    "        new_med = {}\n",
    "        for key,val in medium.items():\n",
    "            if key in ex_ids:\n",
    "                new_med[key] = val\n",
    "        self.model.medium = new_med\n",
    "        self._reset_bounds()\n",
    "        \n",
    "    def optimize(self, disp=False):\n",
    "        sol = linprog(-self.objective_c, A_eq=self.stoichiometry_matrix, b_eq=np.zeros(self.num_metabolites), bounds=self.bounds, method=\"highs\", options={\"disp\":disp})\n",
    "        sol[\"fun\"] = -sol[\"fun\"] # As we have to minimize\n",
    "        return sol \n",
    "    \n",
    "    def slim_optimize(self, disp=False):\n",
    "        sol = self.optimize(disp=disp)\n",
    "        return sol[\"fun\"]\n",
    "\n",
    "    def summary(self):\n",
    "        sol = self.optimize()\n",
    "        flux = sol[\"x\"]\n",
    "        ex_ids = [ex.id for ex in self.exchanges]\n",
    "        fluxes = []\n",
    "        for ex in ex_ids:\n",
    "            idx = self.rec_id[ex]\n",
    "            fluxes.append(flux[idx])\n",
    "        summary_df = pd.DataFrame({\"Exchange reaction\": ex_ids, \"Flux\": fluxes})\n",
    "        summary_df.sort_values([\"Flux\"], inplace=True)\n",
    "        return summary_df\n",
    "\n",
    "    def _reset_bounds(self):\n",
    "        self.bounds = []\n",
    "        for rec in self.model.reactions:\n",
    "            self.bounds.append((rec.lower_bound, rec.upper_bound))\n",
    "\n",
    "    def __add__(self, model2):\n",
    "        \"\"\" Adding another model creates a community model \"\"\"\n",
    "        return CommunityModel([self,model2], [1.,1.])\n",
    "    "
   ]
  },
  {
   "source": [
    "model1 = Model(model_DP, BIOMASS_DP)\n",
    "model2 = Model(model_SA, BIOMASS_SA)"
   ],
   "cell_type": "code",
   "metadata": {},
   "execution_count": 17,
   "outputs": []
  },
  {
   "source": [
    "$$ \\max \\sum_{i=1}^n \\theta_i \\text{ subject to }$$\n",
    "\n",
    "$$ SV = 0, V_{j,min} \\leq V_j \\leq V_{j, max}, V_{Growth} \\geq MBR/10 , V_{i, COOPM} + V_{i, min}\\theta_i \\geq V_{i,\\min} $$"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mip import Model, xsum, maximize, BINARY\n",
    "\n",
    "S1 = model1.stoichiometry_matrix.todense()\n",
    "S1_dict = model1.rec_id \n",
    "bounds1 = model1.bounds \n",
    "obj1 = np.where(model1.objective_c > 0)[0][0]\n",
    "\n",
    "S2 = model2.stoichiometry_matrix.todense()\n",
    "S2_dict = model2.rec_id \n",
    "bounds2 = model2.bounds\n",
    "obj2 = np.where(model2.objective_c > 0)[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "snm3_med = dict(list(model_DP.medium.items()) + list(model_SA.medium.items()))"
   ]
  },
  {
   "source": [
    "We will add following constraints:\n",
    "* All EX_ile__L_e used by DP must be produced by SA.\n",
    "* All EX_met__L_E used by DP must be produced by SA."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mip import Model, xsum, maximize, BINARY\n",
    "\n",
    "comm_model = Model()\n",
    "\n",
    "# Shuttel reactions\n",
    "x_sh = []\n",
    "id1 = []\n",
    "id2 = []\n",
    "x_sh_dict = {}\n",
    "for key, val in snm3_med.items():\n",
    "    x = comm_model.add_var(lb=-val, ub=1000)\n",
    "    x_sh +=[x]\n",
    "    x_sh_dict[key] = x\n",
    "    if key in S1_dict:\n",
    "        id1 += [S1_dict[key]]\n",
    "    else:\n",
    "        id1 += [None]\n",
    "    if key in S2_dict:\n",
    "        id2 += [S2_dict[key]]\n",
    "    else:\n",
    "        id2 += [None]\n",
    "\n",
    "# Flux first model\n",
    "x1 = []\n",
    "for i, (lb, ub) in enumerate(bounds1): \n",
    "    x1 += [comm_model.add_var(lb = lb, ub=ub)]\n",
    "\n",
    "# Flux second model\n",
    "x2 = []\n",
    "for i, (lb, ub) in enumerate(bounds2): \n",
    "    x2 += [comm_model.add_var(lb = lb, ub=ub)]\n",
    "# Stoichiometry\n",
    "for i in range(S1.shape[0]):\n",
    "    comm_model.add_constr(xsum(S1[i,j]*x1[j] for j in range(S1.shape[1])) == 0)\n",
    "\n",
    "for i in range(S2.shape[0]):\n",
    "    comm_model.add_constr(xsum(S2[i,j]*x2[j] for j in range(S2.shape[1])) == 0)\n",
    "\n",
    "ids_ile = (0,0)\n",
    "ids_met = (0,0)\n",
    "# Shuttel constraints\n",
    "for i,key in enumerate(snm3_med):\n",
    "    # if \"EX_ile__L_e\" == key:\n",
    "    #     idx1 = id1[i]\n",
    "    #     idx2 = id2[i]\n",
    "    #     ids_ile = (idx1, idx2)\n",
    "        \n",
    "    #     # Sa must produce it...\n",
    "    #     comm_model.add_constr(x1[idx1] + x2[idx2] == 0)\n",
    "    # elif \"EX_met__L_e\" == key:\n",
    "    #     idx1 = id1[i]\n",
    "    #     idx2 = id2[i]\n",
    "    #     ids_met = (idx1, idx2)\n",
    "    #     # Sa must produce it\n",
    "    #     comm_model.add_constr(x1[idx1] + x2[idx2] == 0)\n",
    "    if id1[i] is not None and id2[i] is not None:\n",
    "        idx1 = id1[i]\n",
    "        idx2 = id2[i]\n",
    "        comm_model.add_constr(-x_sh[i] + x1[idx1] + x2[idx2] == 0)\n",
    "    elif id1[i] is not None:\n",
    "        idx = id1[i]\n",
    "        comm_model.add_constr(-x_sh[i] + x1[idx] == 0)\n",
    "    else:\n",
    "        idx = id2[i]\n",
    "        comm_model.add_constr(-x_sh[i] + x2[idx] == 0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_exchange_flux(S_dict, x):\n",
    "    dic ={}\n",
    "    for key, val in S_dict.items():\n",
    "        if \"EX_\" in key:\n",
    "            dic[key] = x[val].x\n",
    "    return dic\n",
    "\n",
    "def get_medium(x):\n",
    "    dic = {}\n",
    "    for key, val in x.items():\n",
    "        dic[key] = val.lb\n",
    "    return dic"
   ]
  },
  {
   "source": [
    "First looking at the convex combinaiton"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "2.5586946126135386\n0.0\n"
     ]
    }
   ],
   "source": [
    "comm_model.objective = maximize(x1[obj1] + x2[obj2])\n",
    "comm_model.optimize()\n",
    "print(x2[obj2].x)\n",
    "print(x1[obj1].x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "0.23934339136846774\n1.723758259485812\n0.04620056215860694\n"
     ]
    }
   ],
   "source": [
    "# SA dominates\n",
    "w1 = 0.85\n",
    "w2 = 0.15\n",
    "\n",
    "comm_model.objective = maximize(w1*x1[obj1] + w2*x2[obj2])\n",
    "comm_model.optimize()\n",
    "print(x1[obj1].x)\n",
    "print(x2[obj2].x)\n",
    "minMBR = (w1*x1[obj1].x + w2*x2[obj2].x)/10\n",
    "print(minMBR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [],
   "source": [
    "shared_keys = set(min_medium_DP.keys()).intersection(set(min_medium_SA.keys()))\n",
    "shared_medium = dict()\n",
    "for key in shared_keys:\n",
    "    shared_medium[key] = min(min_medium_DP[key], min_medium_SA[key])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "{'EX_his__L_e': 0.002590591762669045,\n",
       " 'EX_fe2_e': 0.0021067642024968536,\n",
       " 'EX_ribflv_e': 1.2195781272038601e-05,\n",
       " 'EX_k_e': 0.0015045698519940388,\n",
       " 'EX_trp__L_e': 0.0015543605265616381,\n",
       " 'EX_cl_e': 0.0012670655721661631,\n",
       " 'EX_glc__D_e': 0.9999999999999942,\n",
       " 'EX_o2_e': 0.012577861281038893,\n",
       " 'EX_mn2_e': 0.00016836210550996272,\n",
       " 'EX_arg__L_e': 0.07486251744579463,\n",
       " 'EX_nac_e': 2.162455164908833e-05,\n",
       " 'EX_leu__L_e': 0.012319680565634147,\n",
       " 'EX_mg2_e': 0.0021116906637898515,\n",
       " 'EX_phe__L_e': 0.008836773254827965,\n",
       " 'EX_so4_e': 0.0006123728011132153,\n",
       " 'EX_lys__L_e': 0.009383669275995604,\n",
       " 'EX_zn2_e': 8.290170544867465e-05,\n",
       " 'EX_na1_e': 0.0025587956488567903,\n",
       " 'EX_thm_e': 6.0978906360193005e-06,\n",
       " 'EX_cys__L_e': 0.0026543116397897693}"
      ]
     },
     "metadata": {},
     "execution_count": 185
    }
   ],
   "source": [
    "shared_medium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [],
   "source": [
    "for id, x in x_sh_dict.items():\n",
    "    if id in shared_medium:\n",
    "        x.lb = -min_medium_DP[id]\n",
    "    else:\n",
    "        x.lb = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "0.0\n0.001286150112519139\n"
     ]
    }
   ],
   "source": [
    "comm_model.objective = maximize(w1*x1[obj1] + w2*x2[obj2])\n",
    "comm_model.optimize()\n",
    "print(x1[obj1].x)\n",
    "print(x2[obj2].x)\n",
    "minMBR = (w1*x1[obj1].x + w2*x2[obj2].x)/10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "EX_glc__D_e -0.868651696924108\nOther exchange EX_glc__D_e -0.1313483030758862\nEX_arg__L_e -0.03287324118981504\nOther exchange EX_arg__L_e -0.0003861381933806529\nEX_o2_e -0.012577861281038894\nOther exchange EX_o2_e 0.0\nEX_glu__L_e -0.031813319336022945\nOther exchange EX_glu__L_e 0.031813319336022945\nEX_ile__L_e -0.025155722562124984\nOther exchange EX_ile__L_e 0.02515572256212499\n"
     ]
    }
   ],
   "source": [
    "interchange = {\"DP_flux\":[], \"SA_flux\":[]}\n",
    "index = list()\n",
    "for key, val in get_exchange_flux(S1_dict,x1).items():\n",
    "    other_echange = get_exchange_flux(S2_dict,x2)\n",
    "    if val < 0:\n",
    "        print(key, val)\n",
    "        index += [key]\n",
    "        if key in other_echange:\n",
    "            print(\"Other exchange\",key, other_echange[key])"
   ]
  },
  {
   "source": [
    "This is expected as the medium does not contain 26dap and SA cannot produce it, thus we will add it to the medium."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "EX_glc__D_e      1.000000\n",
       "EX_leu__L_e      0.012320\n",
       "EX_4abz_e        0.000018\n",
       "EX_cl_e          0.001659\n",
       "EX_pi_e          0.024136\n",
       "EX_ribflv_e      0.000012\n",
       "EX_gly_e         0.028193\n",
       "EX_thr__L_e      0.006937\n",
       "EX_arg__L_e      0.074863\n",
       "EX_lys__L_e      0.009384\n",
       "EX_k_e           0.001505\n",
       "EX_pro__L_e      0.006045\n",
       "EX_ca2_e         0.001468\n",
       "EX_mg2_e         0.002420\n",
       "EX_mn2_e         0.001071\n",
       "EX_cobalt2_e     0.000998\n",
       "EX_zn2_e         0.000899\n",
       "EX_cu2_e         0.000926\n",
       "EX_o2_e          0.012578\n",
       "EX_glu__L_e      0.017145\n",
       "EX_fe2_e         0.002107\n",
       "EX_his__L_e      0.002591\n",
       "EX_ile__L_e      0.007944\n",
       "EX_met__L_e      0.004209\n",
       "EX_so4_e         0.000612\n",
       "EX_trp__L_e      0.001554\n",
       "EX_val__L_e      0.011721\n",
       "EX_thm_e         0.000006\n",
       "EX_cys__L_e      0.002654\n",
       "EX_phe__L_e      0.008837\n",
       "EX_na1_e         0.002559\n",
       "EX_ni2_e         0.001002\n",
       "EX_nac_e         0.000336\n",
       "EX_26dap__M_e    0.002734\n",
       "dtype: float64"
      ]
     },
     "metadata": {},
     "execution_count": 189
    }
   ],
   "source": [
    "min_medium_DP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_medium(medium):\n",
    "    for id, x in x_sh_dict.items():\n",
    "        if id in medium:\n",
    "            x.lb = -min_medium_DP[id]\n",
    "        else:\n",
    "            x.lb = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_medium(shared_medium)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "0.0\n0.001286150112519139\n"
     ]
    }
   ],
   "source": [
    "comm_model.objective = maximize(w1*x1[obj1] + w2*x2[obj2])\n",
    "comm_model.optimize()\n",
    "print(x1[obj1].x)\n",
    "print(x2[obj2].x)\n",
    "minMBR = (w1*x1[obj1].x + w2*x2[obj2].x)/10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import chain, combinations\n",
    "\n",
    "def powerset(iterable):\n",
    "    \"powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)\"\n",
    "    s = list(iterable)\n",
    "    return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {},
   "outputs": [],
   "source": [
    "missing = set(min_medium_DP.keys()).difference(set(shared_keys))\n",
    "missing_combinations = powerset(missing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "('EX_thr__L_e', 'EX_cu2_e', 'EX_ni2_e', 'EX_pro__L_e', 'EX_ca2_e', 'EX_cobalt2_e') 0.00709836307723159\n"
     ]
    }
   ],
   "source": [
    "best_medium=None\n",
    "for combi in missing_combinations:\n",
    "    medium = shared_medium.copy()\n",
    "    medium[\"EX_26dap__M_e\"] = 0.002734\n",
    "    for ex in combi:\n",
    "        medium[ex] = min_medium_DP[ex]\n",
    "    set_medium(medium)\n",
    "    comm_model.objective = maximize(x1[obj1])\n",
    "    comm_model.optimize()\n",
    "    if x1[obj1].x > 0:\n",
    "        if best_medium is None:\n",
    "            print(combi, x1[obj1].x)\n",
    "            best_medium = medium\n",
    "        elif len(medium) < len(best_medium):\n",
    "            print(combi, x1[obj1].x)\n",
    "            best_medium = medium \n",
    "        else:\n",
    "            pass"
   ]
  },
  {
   "source": [
    "As expected these metabolits are **no exchange reactions for SA**, hence they cannot be produced by it. Thus we will have to substitute them into the medium"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "Hence we will add the following increasients to the medium: 'EX_thr__L_e', 'EX_cu2_e', 'EX_ni2_e', 'EX_pro__L_e', 'EX_ca2_e', 'EX_cobalt2_e'"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "{'EX_his__L_e': 0.002590591762669045,\n",
       " 'EX_fe2_e': 0.0021067642024968536,\n",
       " 'EX_ribflv_e': 1.2195781272038601e-05,\n",
       " 'EX_k_e': 0.0015045698519940388,\n",
       " 'EX_trp__L_e': 0.0015543605265616381,\n",
       " 'EX_cl_e': 0.0012670655721661631,\n",
       " 'EX_glc__D_e': 0.9999999999999942,\n",
       " 'EX_o2_e': 0.012577861281038893,\n",
       " 'EX_mn2_e': 0.00016836210550996272,\n",
       " 'EX_arg__L_e': 0.07486251744579463,\n",
       " 'EX_nac_e': 2.162455164908833e-05,\n",
       " 'EX_leu__L_e': 0.012319680565634147,\n",
       " 'EX_mg2_e': 0.0021116906637898515,\n",
       " 'EX_phe__L_e': 0.008836773254827965,\n",
       " 'EX_so4_e': 0.0006123728011132153,\n",
       " 'EX_lys__L_e': 0.009383669275995604,\n",
       " 'EX_zn2_e': 8.290170544867465e-05,\n",
       " 'EX_na1_e': 0.0025587956488567903,\n",
       " 'EX_thm_e': 6.0978906360193005e-06,\n",
       " 'EX_cys__L_e': 0.0026543116397897693,\n",
       " 'EX_26dap__M_e': 0.002734,\n",
       " 'EX_thr__L_e': 0.0069370205460978955,\n",
       " 'EX_cu2_e': 0.0009257250408242616,\n",
       " 'EX_ni2_e': 0.0010022613010017905,\n",
       " 'EX_pro__L_e': 0.006044704997960757,\n",
       " 'EX_ca2_e': 0.0014677908938624288,\n",
       " 'EX_cobalt2_e': 0.0009981830860061644}"
      ]
     },
     "metadata": {},
     "execution_count": 205
    }
   ],
   "source": [
    "best_medium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_medium(best_medium)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "0.007098367882456253\n0.0009628255517842001\n"
     ]
    }
   ],
   "source": [
    "comm_model.objective = maximize(w1*x1[obj1] + w2*x2[obj2])\n",
    "comm_model.optimize()\n",
    "print(x1[obj1].x)\n",
    "print(x2[obj2].x)\n",
    "minMBR = (w1*x1[obj1].x + w2*x2[obj2].x)/10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "EX_h2o_e -0.029945258591327507\nOther exchange EX_h2o_e 0.4904757172624663\nEX_arg__L_e -0.025080738842006508\nOther exchange EX_arg__L_e -0.049781778603788115\nEX_his__L_e -5.800741764403502e-05\nOther exchange EX_his__L_e -0.0006512474209312142\nEX_k_e -0.00017878707671080808\nOther exchange EX_k_e -0.0003782329774385043\nEX_leu__L_e -0.0003023377080732263\nOther exchange EX_leu__L_e -0.0030970376385353383\nEX_lys__L_e -0.0003470738841987999\nOther exchange EX_lys__L_e -0.0023589553950281437\nEX_cl_e -4.767912132435359e-06\nOther exchange EX_cl_e -0.00041712313264840407\nEX_cys__L_e -0.001218179949249762\nOther exchange EX_cys__L_e -0.0014361316905400073\nEX_mg2_e -7.946199278875004e-06\nOther exchange EX_mg2_e -0.0006084454401063097\nEX_mn2_e -6.335392130740038e-07\nOther exchange EX_mn2_e -0.0002691807716321314\nEX_mobd_e -6.7397788630923146e-09\nEX_na1_e -0.0019155414352746666\nOther exchange EX_na1_e -0.0006432542135821236\nEX_fe3_e -1.4209379493231231e-05\nEX_nac_e -8.137231410366602e-08\nOther exchange EX_nac_e -8.453604246182182e-05\nEX_orn_e -0.04774844691914762\nOther exchange EX_orn_e 0.04774844691914762\nEX_phe__L_e -0.00011004099528810623\nOther exchange EX_phe__L_e -0.0022214715087457222\nEX_glc__D_e -0.1671301805911462\nOther exchange EX_glc__D_e -0.8328698194088482\nEX_ribflv_e -2.874088832447543e-07\nOther exchange EX_ribflv_e -3.0658906041619967e-06\nEX_so4_e -3.972618226661611e-06\nOther exchange EX_so4_e -0.0001539440575393945\nEX_trp__L_e -1.4769673058165723e-05\nOther exchange EX_trp__L_e -0.0003907498273975913\nEX_tyr__L_e -7.345820376565916e-05\nOther exchange EX_tyr__L_e 0.0\nEX_thm_e -4.294201961341676e-07\nOther exchange EX_thm_e -1.5329453320490266e-06\nEX_thym_e -3.3100491917596106e-05\nEX_zn2_e -3.1195547880598716e-07\nOther exchange EX_zn2_e -0.0002260891684903279\nEX_dcaACP_e -7.803533702377229e-06\n"
     ]
    }
   ],
   "source": [
    "interchange = {\"DP_flux\":[], \"SA_flux\":[]}\n",
    "index = list()\n",
    "for key, val in get_exchange_flux(S2_dict,x2).items():\n",
    "    other_echange = get_exchange_flux(S1_dict,x1)\n",
    "    if val < 0:\n",
    "        print(key, val)\n",
    "        index += [key]\n",
    "        if key in other_echange:\n",
    "            print(\"Other exchange\",key, other_echange[key])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}