community_model_shuttle_Exchange_14.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 245,
   "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 json"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": "code",
   "execution_count": 246,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Growth:  0.2823653925322476\n",
      "Growth:  2.558694612613393\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": 247,
   "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": 248,
   "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>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>"
      ],
      "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"
      ]
     },
     "execution_count": 248,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "snm3 = pd.read_csv(\"SNM3.csv\", sep =\";\")\n",
    "snm3.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {},
   "outputs": [],
   "source": [
    "BIOMASS_DP = \"Growth\" \n",
    "BIOMASS_SA = \"BIOMASS_iYS_wild_type\"\n",
    "models = [model_DP.copy(), model_SA.copy()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Defining exchange media"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of united exchange reactions:  111\n"
     ]
    }
   ],
   "source": [
    "ex_DP_ids = set(map(lambda x:x.id, model_DP.exchanges))\n",
    "ex_SA_ids = set(map(lambda x:x.id, model_SA.exchanges))\n",
    "ex_union_ids = ex_DP_ids.union(ex_SA_ids)\n",
    "print(\"Number of united exchange reactions: \", len(ex_union_ids))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "metadata": {},
   "outputs": [],
   "source": [
    "fva_DP = cobra.flux_analysis.flux_variability_analysis(model_DP)\n",
    "fva_SA = cobra.flux_analysis.flux_variability_analysis(model_SA)                  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining medium with all 111 total exchange reaction\n",
    "#def compute_exchange(fva_DP, fva_SA):\n",
    "#    media_exchange_DP = dict()\n",
    "#    media_exchange_SA = dict()\n",
    "#    for ex in ex_union_ids:\n",
    "#        flux = - 10\n",
    "#        if flux < 0:\n",
    "#            media_exchange_DP[ex] = abs(flux)\n",
    "#            media_exchange_SA[ex] = abs(flux)\n",
    "#        \n",
    "#    return media_exchange_DP, media_exchange_SA  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_exchange(fva_DP, fva_SA):\n",
    "    media_exchange_DP = dict()\n",
    "    media_exchange_SA = dict()\n",
    "    for ex in ex_union_ids:\n",
    "        if ex in fva_DP.index and ex in fva_SA.index:\n",
    "            if not (fva_DP.loc[ex].minimum ==0 and fva_DP.loc[ex].maximum == 0):\n",
    "                if not (fva_SA.loc[ex].minimum ==0 and fva_SA.loc[ex].maximum == 0):\n",
    "                    flux = - 10\n",
    "                    if flux < 0:\n",
    "                        media_exchange_DP[ex] = abs(flux)\n",
    "                        media_exchange_SA[ex] = abs(flux)\n",
    "        elif ex in fva_DP.index:\n",
    "            # Take just the minimum if not common\n",
    "            flux = fva_DP.loc[ex].minimum\n",
    "            if flux < 0:\n",
    "                media_exchange_DP[ex] = abs(flux)\n",
    "        elif ex in fva_SA.index:\n",
    "            # Take just the minimum if not common\n",
    "            flux = fva_SA.loc[ex].minimum\n",
    "            if flux < 0:\n",
    "                media_exchange_SA[ex] = abs(flux)\n",
    "        else:\n",
    "            print(\"error\")\n",
    "    return media_exchange_DP, media_exchange_SA  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [],
   "source": [
    "media_exchange_DP, media_exchange_SA = compute_exchange(fva_DP,fva_SA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_dict(data, name):\n",
    "    with open(name, 'w' ) as file:\n",
    "        json.dump( data, file )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of media:  44\n",
      "{'EX_ala__L_e': 10.0, 'EX_na1_e': 10, 'EX_lys__L_e': 10, 'EX_phe__L_e': 10, 'EX_pro__L_e': 0.06044704997960771, 'EX_glu__L_e': 10, 'EX_tyr__L_e': 10, 'EX_h_e': 10, 'EX_leu__L_e': 10, 'EX_glyclt_e': 10, 'EX_cu2_e': 0.009257250408242637, 'EX_fe2_e': 10, 'EX_o2_e': 10, 'EX_ile__L_e': 10, 'EX_co2_e': 10, 'EX_arg__L_e': 10, 'EX_cys__L_e': 10, 'EX_thm_e': 10, 'EX_so4_e': 10, 'EX_orn_e': 10, 'EX_nh4_e': 10, 'EX_4abz_e': 0.00018293671908059108, 'EX_ribflv_e': 10, 'EX_gly_e': 10.0, 'EX_mg2_e': 10, 'EX_zn2_e': 10, 'EX_glc__D_e': 10, 'EX_pi_e': 10, 'EX_trp__L_e': 10, 'EX_ca2_e': 0.01467790893862432, 'EX_nac_e': 10, 'EX_fum_e': 10.0, 'EX_thr__L_e': 0.06937020546135564, 'EX_k_e': 10, 'EX_mn2_e': 10, 'EX_val__L_e': 10, 'EX_ni2_e': 0.010022613010017837, 'EX_his__L_e': 10, 'EX_met__L_e': 10, 'EX_h2o_e': 10, 'EX_cobalt2_e': 0.009981830860061667, 'EX_acald_e': 10, 'EX_26dap__M_e': 0.027344801057605284, 'EX_cl_e': 10}\n"
     ]
    }
   ],
   "source": [
    "print(\"Length of media: \", len(media_exchange_DP))\n",
    "print(media_exchange_DP)\n",
    "save_dict(media_exchange_DP, \"exchange_DP.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of media:  37\n",
      "{'EX_na1_e': 10, 'EX_lys__L_e': 10, 'EX_phe__L_e': 10, 'EX_glu__L_e': 10, 'EX_tyr__L_e': 10, 'EX_h_e': 10, 'EX_leu__L_e': 10, 'EX_pyr_e': 10.0, 'EX_glyclt_e': 10, 'EX_fe2_e': 10, 'EX_o2_e': 10, 'EX_urea_e': 10.0, 'EX_ile__L_e': 10, 'EX_co2_e': 10, 'EX_arg__L_e': 10, 'EX_cys__L_e': 10, 'EX_thm_e': 10, 'EX_so4_e': 10, 'EX_orn_e': 10, 'EX_nh4_e': 10, 'EX_ribflv_e': 10, 'EX_ser__L_e': 10.0, 'EX_mg2_e': 10, 'EX_zn2_e': 10, 'EX_glc__D_e': 10, 'EX_pi_e': 10, 'EX_trp__L_e': 10, 'EX_nac_e': 10, 'EX_k_e': 10, 'EX_mn2_e': 10, 'EX_val__L_e': 10, 'EX_his__L_e': 10, 'EX_met__L_e': 10, 'EX_h2o_e': 10, 'EX_acald_e': 10, 'EX_cl_e': 10, 'EX_mobd_e': 1.791086228829375e-05}\n"
     ]
    }
   ],
   "source": [
    "print(\"Length of media: \", len(media_exchange_SA))\n",
    "print(media_exchange_SA)\n",
    "save_dict(media_exchange_SA, \"exchange_SA.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "exchange_SA = json.loads(open(\"exchange_SA.json\").read())\n",
    "exchange_DP = json.loads(open(\"exchange_DP.json\").read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [],
   "source": [
    "exchange = dict()\n",
    "for key, val in exchange_SA.items():\n",
    "    exchange[key] = val \n",
    "for key, val in exchange_DP.items():\n",
    "    exchange[key] = val "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Constructing of community model explicitely"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "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": 261,
   "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",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CommunityModel(Model):\n",
    "    def __init__(self, models, weight):\n",
    "        self.models = models\n",
    "        self.id = \"|\".join([model.id for model in models])\n",
    "        self.shared_exchanges = []\n",
    "        self.weight = weight\n",
    "        for model in models:\n",
    "            for rec in model.exchanges:\n",
    "                if rec.id not in self.shared_exchanges:\n",
    "                    self.shared_exchanges.append(rec.id)\n",
    "        # Create community stoichimetry matrix with shuttel reactions!\n",
    "        self._shifts = [len(self.shared_exchanges)]\n",
    "        for i,model in enumerate(models):\n",
    "            self._shifts.append(self._shifts[i] + model.num_reactions)\n",
    "        S_EX = -np.eye(self._shifts[0])\n",
    "        matrices = [S_EX] + [m.stoichiometry_matrix.todense() for m in models]\n",
    "        S = scipy.linalg.block_diag(*matrices)\n",
    "        self.num_reactions = S.shape[1]\n",
    "        self.num_metabolites = S.shape[0]\n",
    "        for i, id in enumerate(self.shared_exchanges):\n",
    "            for j,model in enumerate(models):\n",
    "                if id in model.rec_id:\n",
    "                    S[i,self._shifts[j] + model.rec_id[id]] = 1\n",
    "        self.stoichiometry_matrix = scipy.sparse.csr_matrix(S)\n",
    "        # Cretae objective:\n",
    "        self._weighted_objective(weight)\n",
    "        # Create bounds\n",
    "        self._reset_bounds()\n",
    "    \n",
    "    @property\n",
    "    def reactions(self):\n",
    "        return [model.reactions for model in self.models]\n",
    "    @property\n",
    "    def exchanges(self):\n",
    "        return self.shared_exchanges\n",
    "    @property\n",
    "    def metabolites(self):\n",
    "        return [model.metabolites for model in self.models]\n",
    "    @property\n",
    "    def medium(self):\n",
    "        medium = {}\n",
    "        for model in self.models:\n",
    "            for key,val in model.medium:\n",
    "                medium[key] = val\n",
    "        return medium\n",
    "\n",
    "    def _weighted_objective(self, weight):\n",
    "        self.weight = weight\n",
    "        self.objective_c = np.zeros(self._shifts[0])\n",
    "        for i,model in enumerate(self.models):\n",
    "            self.objective_c = np.append(self.objective_c, weight[i]*model.objective_c)\n",
    "    def _reset_bounds(self):\n",
    "        self.bounds = []\n",
    "        for id in self.shared_exchanges:\n",
    "            min_lower_bound = 0\n",
    "            for model in self.models:\n",
    "                if id in model.rec_id:\n",
    "                    rec = model.reactions.get_by_id(id)\n",
    "                    if rec.lower_bound < min_lower_bound:\n",
    "                        min_lower_bound = rec.lower_bound \n",
    "            self.bounds.append((min_lower_bound, 1000))\n",
    "        for model in self.models:\n",
    "            self.bounds += model.bounds\n",
    "\n",
    "    def set_medium(self, medium):\n",
    "        for model in self.models:\n",
    "            model.set_medium(medium)\n",
    "        self._reset_bounds()\n",
    "\n",
    "    def summary(self):\n",
    "        sol = self.optimize()\n",
    "        flux = sol[\"x\"]\n",
    "        ex_ids = self.shared_exchanges\n",
    "        ex_flux = flux[:len(ex_ids)]\n",
    "        df_ex = pd.DataFrame({\"Exchanges\": ex_ids, \"Flux\": ex_flux})\n",
    "        df_ex.sort_values([\"Flux\"], inplace=True)\n",
    "        for i,model in enumerate(self.models):\n",
    "            shuttel_ids = [ex.id for ex in model.exchanges]\n",
    "            id = str(model.id) + \" Shuttel Flux\"\n",
    "            df_ex[id] = 0.\n",
    "            for sh in shuttel_ids:\n",
    "                idx = model.rec_id[sh]\n",
    "                df_ex[id][ex_ids.index(sh)] = flux[self._shifts[i] +idx]\n",
    "        \n",
    "        return df_ex\n",
    "    \n",
    "    def set_weights(self, weight):\n",
    "        self._weighted_objective(weight)\n",
    "    \n",
    "    def get_model_growths(self):\n",
    "        mask = self.objective_c != 0\n",
    "        sol = self.optimize()\n",
    "        flux = sol[\"x\"]\n",
    "        return flux[mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 263,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Growth DP:  0.28236539253224774\n",
      "Growth SA:  2.5586946126118373\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Sparse constraint matrix detected; setting 'sparse':True.\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",
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Unknown solver options: sparse\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"
     ]
    }
   ],
   "source": [
    "model_DP = Model(model_DP, BIOMASS_DP)\n",
    "print(\"Growth DP: \",model_DP.slim_optimize())\n",
    "model_SA = Model(model_SA, BIOMASS_SA)\n",
    "print(\"Growth SA: \",model_SA.slim_optimize())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weights 1:1:  2.5586946126136296\n",
      "SA growth:  2.5586946126136296\n",
      "DP growth:  -0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Sparse constraint matrix detected; setting 'sparse':True.\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",
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Unknown solver options: sparse\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"
     ]
    }
   ],
   "source": [
    "# Community model\n",
    "model = model_SA + model_DP\n",
    "print(\"Weights 1:1: \", model.slim_optimize())\n",
    "growths = model.get_model_growths()\n",
    "print(\"SA growth: \", growths[0])\n",
    "print(\"DP growth: \", growths[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.set_medium(exchange)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Weights 1:1:  2.2982805284808725\n",
      "SA growth:  2.2982805284808725\n",
      "DP growth:  -0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Sparse constraint matrix detected; setting 'sparse':True.\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",
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Unknown solver options: sparse\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"
     ]
    }
   ],
   "source": [
    "model = model_SA + model_DP\n",
    "print(\"Weights 1:1: \", model.slim_optimize())\n",
    "growths = model.get_model_growths()\n",
    "print(\"SA growth: \", growths[0])\n",
    "print(\"DP growth: \", growths[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Sparse constraint matrix detected; setting 'sparse':True.\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",
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Unknown solver options: sparse\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",
      "<ipython-input-262-65ea33eaf871>:83: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_ex[id][ex_ids.index(sh)] = flux[self._shifts[i] +idx]\n",
      "<ipython-input-262-65ea33eaf871>:83: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_ex[id][ex_ids.index(sh)] = flux[self._shifts[i] +idx]\n"
     ]
    },
    {
     "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>Exchanges</th>\n",
       "      <th>Flux</th>\n",
       "      <th>iYS854 Shuttel Flux</th>\n",
       "      <th>DP_83VPs_KB5 Shuttel Flux</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>EX_pyr_e</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62</th>\n",
       "      <td>EX_urea_e</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>EX_h_e</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>51</th>\n",
       "      <td>EX_glc__D_e</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>EX_acald_e</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>98</th>\n",
       "      <td>EX_ser__D_e</td>\n",
       "      <td>10.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>EX_etoh_e</td>\n",
       "      <td>11.509173</td>\n",
       "      <td>11.509173</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>EX_lac__L_e</td>\n",
       "      <td>22.948626</td>\n",
       "      <td>22.948626</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>EX_h2o_e</td>\n",
       "      <td>23.016521</td>\n",
       "      <td>33.016521</td>\n",
       "      <td>-10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>EX_co2_e</td>\n",
       "      <td>25.259855</td>\n",
       "      <td>25.259855</td>\n",
       "      <td>-0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>111 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      Exchanges       Flux  iYS854 Shuttel Flux  DP_83VPs_KB5 Shuttel Flux\n",
       "52     EX_pyr_e -10.000000           -10.000000                        0.0\n",
       "62    EX_urea_e -10.000000           -10.000000                        0.0\n",
       "12       EX_h_e -10.000000           -10.000000                       -0.0\n",
       "51  EX_glc__D_e -10.000000           -10.000000                       -0.0\n",
       "5    EX_acald_e -10.000000           -10.000000                       -0.0\n",
       "..          ...        ...                  ...                        ...\n",
       "98  EX_ser__D_e  10.000000             0.000000                       10.0\n",
       "32    EX_etoh_e  11.509173            11.509173                        0.0\n",
       "19  EX_lac__L_e  22.948626            22.948626                        0.0\n",
       "10     EX_h2o_e  23.016521            33.016521                      -10.0\n",
       "23     EX_co2_e  25.259855            25.259855                       -0.0\n",
       "\n",
       "[111 rows x 4 columns]"
      ]
     },
     "execution_count": 267,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = model.summary()\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv(\"Summary_exchange_shuttels_Exchange_weight_1_1_RM.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 269,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "weights 1:10  4.032449748658346\n",
      "SA growth:  1.2087958233205274\n",
      "DP growth:  0.28236539253378184\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Sparse constraint matrix detected; setting 'sparse':True.\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",
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Unknown solver options: sparse\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"
     ]
    }
   ],
   "source": [
    "model.set_weights([1,10])\n",
    "print(\"weights 1:10 \", model.slim_optimize())\n",
    "growths = model.get_model_growths()\n",
    "print(\"SA growth: \", growths[0])\n",
    "print(\"DP growth: \", growths[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Sparse constraint matrix detected; setting 'sparse':True.\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",
      "<ipython-input-261-5ceee4aecab6>:49: OptimizeWarning: Unknown solver options: sparse\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",
      "<ipython-input-262-65ea33eaf871>:83: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_ex[id][ex_ids.index(sh)] = flux[self._shifts[i] +idx]\n",
      "<ipython-input-262-65ea33eaf871>:83: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_ex[id][ex_ids.index(sh)] = flux[self._shifts[i] +idx]\n"
     ]
    },
    {
     "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>Exchanges</th>\n",
       "      <th>Flux</th>\n",
       "      <th>iYS854 Shuttel Flux</th>\n",
       "      <th>DP_83VPs_KB5 Shuttel Flux</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>EX_o2_e</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-9.911929e+00</td>\n",
       "      <td>-0.088071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>EX_cys__L_e</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-1.200710e-01</td>\n",
       "      <td>-9.879929</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>EX_h_e</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-1.000000e+01</td>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>51</th>\n",
       "      <td>EX_glc__D_e</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>5.314926e-11</td>\n",
       "      <td>-10.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>56</th>\n",
       "      <td>EX_ser__L_e</td>\n",
       "      <td>-10.000000</td>\n",
       "      <td>-1.000000e+01</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>103</th>\n",
       "      <td>EX_h2s_e</td>\n",
       "      <td>9.853386</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>9.853386</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>EX_etoh_e</td>\n",
       "      <td>10.775537</td>\n",
       "      <td>1.077554e+01</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>EX_lac__L_e</td>\n",
       "      <td>12.730453</td>\n",
       "      <td>1.273045e+01</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>EX_co2_e</td>\n",
       "      <td>26.968431</td>\n",
       "      <td>2.503038e+01</td>\n",
       "      <td>1.938053</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>98</th>\n",
       "      <td>EX_ser__D_e</td>\n",
       "      <td>28.302800</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>28.302800</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>111 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       Exchanges       Flux  iYS854 Shuttel Flux  DP_83VPs_KB5 Shuttel Flux\n",
       "42       EX_o2_e -10.000000        -9.911929e+00                  -0.088071\n",
       "24   EX_cys__L_e -10.000000        -1.200710e-01                  -9.879929\n",
       "12        EX_h_e -10.000000        -1.000000e+01                  -0.000000\n",
       "51   EX_glc__D_e -10.000000         5.314926e-11                 -10.000000\n",
       "56   EX_ser__L_e -10.000000        -1.000000e+01                   0.000000\n",
       "..           ...        ...                  ...                        ...\n",
       "103     EX_h2s_e   9.853386         0.000000e+00                   9.853386\n",
       "32     EX_etoh_e  10.775537         1.077554e+01                   0.000000\n",
       "19   EX_lac__L_e  12.730453         1.273045e+01                   0.000000\n",
       "23      EX_co2_e  26.968431         2.503038e+01                   1.938053\n",
       "98   EX_ser__D_e  28.302800         0.000000e+00                  28.302800\n",
       "\n",
       "[111 rows x 4 columns]"
      ]
     },
     "execution_count": 270,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = model.summary()\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv(\"Summary_exchange_shuttels_Exchange_weight_10_1_RM.csv\")"
   ]
  }
 ],
 "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.8.5"
  },
  "metadata": {
   "interpreter": {
    "hash": "a8f61be024eba58adef938c9aa1e29e02cb3dece83a5348b1a2dafd16a070453"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}