Generate SOY tissue specific models 3-backup.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modelling metabolism in SOY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating SOY core model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from cobra import io\n",
    "from cobra.core import Metabolite, Reaction\n",
    "from cobra import flux_analysis\n",
    "from Functions import *\n",
    "from IPython import display\n",
    "import logging\n",
    "logging.basicConfig()\n",
    "logger = logging.getLogger('logger')\n",
    "\n",
    "#import sbml file\n",
    "fname = \"C:\\\\Users\\\\FFAR\\\\Downloads\\\\plantcoremetabolism-model\\\\PlantCoreMetabolism_v1_3_2.xml\"\n",
    "model = io.read_sbml_model(fname)\n",
    "model.solver='glpk'\n",
    "display.clear_output()\n",
    "\n",
    "from sweetlovegroup import transform\n",
    "model = transform.fixModelCompatibilityIssueCobra015(model,fname)\n",
    "\n",
    "\n",
    "#change Biomass_tx to Protein_biomass\n",
    "rxn = model.reactions.get_by_id(\"Biomass_tx\")\n",
    "rxn.id = \"Protein_biomass_tx\"\n",
    "for met in [\"K_b\",\"Ca_b\",\"Mg_b\"]:\n",
    "    met  = model.metabolites.get_by_id(met)\n",
    "    coeff = rxn.metabolites.get(met)\n",
    "    rxn.add_metabolites({met:-1*coeff})\n",
    "met = Metabolite(\"Protein_b\",name=\"Protein_biomass\")\n",
    "formula_dict = rxn.check_mass_balance()\n",
    "met.formula = \"\".join([atom+str(formula_dict[atom]*-1) for atom in formula_dict.keys() if atom != \"charge\"])\n",
    "met.charge = formula_dict[\"charge\"]*-1\n",
    "met.compartment = \"b\"\n",
    "rxn.add_metabolites({met:1})\n",
    "\n",
    "\n",
    "for rxn in model.reactions:\n",
    "    if rxn.lower_bound == -1000:\n",
    "        rxn.lower_bound = -3000\n",
    "    if rxn.upper_bound == 1000:\n",
    "        rxn.upper_bound = 3000\n",
    "\n",
    "\n",
    "model.metabolites.get_by_id(\"aL_1_PHOSPHATIDYL_GLYCEROL_P_m\").remove_from_model()\n",
    "model.metabolites.get_by_id(\"aL_1_PHOSPHATIDYL_GLYCEROL_P_p\").remove_from_model()\n",
    "\n",
    "model.reactions.get_by_id(\"PGPPHOSPHA_RXN_m\").add_metabolites({model.metabolites.get_by_id(\"L_1_PHOSPHATIDYL_GLYCEROL_P_m\"):-0.03,\n",
    "                                                               model.metabolites.get_by_id(\"PROTON_m\"):-0.03})\n",
    "model.reactions.get_by_id(\"PHOSPHAGLYPSYN_RXN_m\").add_metabolites({model.metabolites.get_by_id(\"L_1_PHOSPHATIDYL_GLYCEROL_P_m\"):0.03,\n",
    "                                                                   model.metabolites.get_by_id(\"PROTON_m\"):0.03})\n",
    "model.reactions.get_by_id(\"PGPPHOSPHA_RXN_p\").add_metabolites({model.metabolites.get_by_id(\"L_1_PHOSPHATIDYL_GLYCEROL_P_p\"):-0.03,\n",
    "                                                               model.metabolites.get_by_id(\"PROTON_p\"):-0.03})\n",
    "model.reactions.get_by_id(\"PHOSPHAGLYPSYN_RXN_p\").add_metabolites({model.metabolites.get_by_id(\"L_1_PHOSPHATIDYL_GLYCEROL_P_p\"):0.03,\n",
    "                                                                   model.metabolites.get_by_id(\"PROTON_p\"):0.03})\n",
    "model.reactions.get_by_id(\"LPG_biosynthesis_c\").add_metabolites({model.metabolites.get_by_id(\"L_1_PHOSPHATIDYL_GLYCEROL_P_p\"):-0.03})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_csv(\"Data/biomass_soy.csv\")\n",
    "\n",
    "FA=[\"PALMITATE_p\",\"CPD_9245_p\",\"CPD_17412_p\",\"CPD_17291_p\",\"STEARIC_ACID_p\",\"OLEATE_CPD_p\",\n",
    "    \"Octadecadienoate_p\",\"LINOLENIC_ACID_p\",\"ARACHIDIC_ACID_p\",\"CPD_16709_p\",\"DOCOSANOATE_p\"]\n",
    "FACP = {\"PALMITATE_p\":\"Palmitoyl_ACPs_p\",\n",
    "        \"CPD_9245_p\":\"Palmitoleoyl_ACP_p\",\n",
    "        \"CPD_17412_p\":\"hexadecadienoate_ACP_p\",\n",
    "        \"CPD_17291_p\":\"hexadecatrienoate_ACP_p\",\n",
    "        \"STEARIC_ACID_p\":\"Stearoyl_ACPs_p\",\n",
    "        \"OLEATE_CPD_p\":\"Oleoyl_ACPs_p\",\n",
    "        \"Octadecadienoate_p\":\"Octadecadienoyl_ACP_p\",\n",
    "        \"LINOLENIC_ACID_p\":\"Octadecatrienoyl_ACP_p\",\n",
    "        \"ARACHIDIC_ACID_p\":\"Arachidoyl_ACPs_p\",\n",
    "        \"CPD_16709_p\":\"Eicosenoyl_ACP_p\",\n",
    "        \"DOCOSANOATE_p\":\"Behenoyl_ACPs_p\"}\n",
    "\n",
    "\n",
    "PLs = [\"ACYL_SN_GLYCEROL_3P_p\",\n",
    "       \"L_PHOSPHATIDATE_p\",\"L_PHOSPHATIDATE_m\",\"DIACYLGLYCEROL_p\",\n",
    "       \"DIACYLGLYCEROL_r\",\"Triacylglycerols_p\",\"PHOSPHATIDYL_CHOLINE_r\",\n",
    "       \"L_1_PHOSPHATIDYL_ETHANOLAMINE_r\",\"L_1_PHOSPHATIDYL_GLYCEROL_p\",\n",
    "       \"L_1_PHOSPHATIDYL_GLYCEROL_P_p\",\"L_1_PHOSPHATIDYL_GLYCEROL_P_m\",\n",
    "       \"L_1_PHOSPHATIDYL_GLYCEROL_m\",\"2_Lysophosphatidylcholines_r\",\n",
    "       \"Lysophosphatidylglycerols_r\",\"CDPDIACYLGLYCEROL_p\",\"CDPDIACYLGLYCEROL_m\",\n",
    "       \"D_Galactosyl_12_diacyl_glycerols_p\",\"Galactosyl_galactosyl_diacyl_glycerols_p\"]\n",
    "\n",
    "\n",
    "for met in PLs:\n",
    "    met=model.metabolites.get_by_id(met)\n",
    "    met.formula=\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use biomass composition to generate biomass equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generateMissingFormula(model,debug=False):\n",
    "    loop_counter = 0\n",
    "    former = 0\n",
    "    for met in model.metabolites:\n",
    "        if met.formula == \"\" or met.formula == \"NA\":\n",
    "            former = former +1\n",
    "    latter = 1\n",
    "    while True:\n",
    "        loop_counter = loop_counter+1\n",
    "        if debug:\n",
    "            print(\"Loop = \"+str(loop_counter))\n",
    "        former = latter\n",
    "        for rxn in model.reactions:\n",
    "            count = 0\n",
    "            for met in rxn.metabolites:\n",
    "                if met.formula==\"\" or met.formula==\"NA\" or met.formula == None:\n",
    "                    if met.formula == \"NA\" or met.formula == None:\n",
    "                        met.formula = \"\"\n",
    "                    count = count + 1\n",
    "                    coeff = rxn.metabolites[met]\n",
    "            if count == 1:\n",
    "                unb = rxn.check_mass_balance()\n",
    "                eqn = rxn.reaction\n",
    "                eqn = \" \"+eqn+\" \"\n",
    "                for met in rxn.metabolites.keys():\n",
    "                    formula = met.formula\n",
    "                    if formula == None:\n",
    "                        formula = \"0\"\n",
    "                        NF_list.add(rxn.id)\n",
    "                    eqn=eqn.replace(\" \"+met.id+\" \",\"(\"+formula+\")\")\n",
    "                if debug:\n",
    "                    print(eqn)\n",
    "                    print(unb)\n",
    "                for met in rxn.metabolites:\n",
    "                    if met.formula == \"\":\n",
    "                        tempForm = \"\"\n",
    "                        for a in sorted(unb.keys()):\n",
    "                            if a==\"charge\" or round(unb[a],2)==0:\n",
    "                                continue\n",
    "                            num = float(abs(unb[a]))/abs(coeff)\n",
    "                            if str(round(num))==str(num):\n",
    "                                tempForm = tempForm+a+str(int(round(num)))\n",
    "                            else:\n",
    "                                tempForm = tempForm+a+str(num)\n",
    "                                #print(a)\n",
    "                                #print(round(num)==num)\n",
    "                                #print(round(num))\n",
    "                                #print(num)\n",
    "                                #print(type(round(num)))\n",
    "                                #print(type(num))\n",
    "                        met.formula = tempForm\n",
    "                        if debug:\n",
    "                            print(met.id)\n",
    "                            print(tempForm)\n",
    "        latter = 0\n",
    "        for met in model.metabolites:\n",
    "            if met.formula == \"\" or met.formula == \"NA\":\n",
    "                latter = latter +1\n",
    "        if former == latter:\n",
    "            break\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run stem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GLN_ph + PROTON_e --> GLN_c + PROTON_c\n",
      "ASN_ph + PROTON_e --> ASN_c + PROTON_c\n",
      "ILE_ph + PROTON_e --> ILE_c + PROTON_c\n",
      "L_ALPHA_ALANINE_ph + PROTON_e --> L_ALPHA_ALANINE_c + PROTON_c\n",
      "GLT_ph + PROTON_e --> GLT_c + PROTON_c\n",
      "PROTON_e + VAL_ph --> PROTON_c + VAL_c\n",
      "GLY_ph + PROTON_e --> GLY_c + PROTON_c\n",
      "PROTON_e + THR_ph --> PROTON_c + THR_c\n",
      "4_AMINO_BUTYRATE_ph + PROTON_e --> 4_AMINO_BUTYRATE_c + PROTON_c\n",
      "LYS_ph + PROTON_e --> LYS_c + PROTON_c\n",
      "PROTON_e + SER_ph --> PROTON_c + SER_c\n",
      "LEU_ph + PROTON_e --> LEU_c + PROTON_c\n",
      "PROTON_e + sSUCROSE_ph --> PROTON_c + sSUCROSE_b\n",
      "L_ASPARTATE_ph + PROTON_e --> L_ASPARTATE_c + PROTON_c\n",
      "GLC_ph + PROTON_e --> GLC_c + PROTON_c\n",
      "FRU_ph + PROTON_e --> FRU_c + PROTON_c\n",
      "PHE_ph + PROTON_e --> PHE_c + PROTON_c\n",
      "ARG_ph + PROTON_e --> ARG_c + PROTON_c\n",
      "HIS_ph + PROTON_e --> HIS_c + PROTON_c\n",
      "CYS_ph + PROTON_e --> CYS_c + PROTON_c\n",
      "PROTON_e + PRO_ph --> PROTON_c + PRO_c\n",
      "MET_ph + PROTON_e --> MET_c + PROTON_c\n",
      "PROTON_e + TYR_ph --> PROTON_c + TYR_c\n",
      "PROTON_e + TRP_ph --> PROTON_c + TRP_c\n"
     ]
    }
   ],
   "source": [
    "temp  = model.copy()\n",
    "\n",
    "for met in temp.reactions.Phloem_output_tx.metabolites.keys():\n",
    "    met2 = met.copy()\n",
    "    if met.id==\"sSUCROSE_b\":\n",
    "        met2.id = met.id.replace(\"_b\",\"_ph\")\n",
    "    elif \"PROTON\" in met.id:\n",
    "        continue\n",
    "    else:\n",
    "        met2.id = met.id.replace(\"_c\",\"_ph\")\n",
    "    met2.compartment = \"ph\"\n",
    "    temp.add_metabolites(met2)\n",
    "    \n",
    "    rxn = Reaction(met2.id.replace(\"_ph\",\"_phloem_uptake\"),name=met2.id.replace(\"_ph\",\"_phloem_uptake\"))\n",
    "    rxn.add_metabolites({met2:-1,temp.metabolites.get_by_id(\"PROTON_e\"):-1,met:1,temp.metabolites.get_by_id(\"PROTON_c\"):1})\n",
    "    rxn.lower_bound = 0\n",
    "    rxn.upper_bound = 1000\n",
    "    temp.add_reaction(rxn)\n",
    "    \n",
    "    print(rxn.reaction)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generateStemModel(model):\n",
    "    for met in model.reactions.Phloem_output_tx.metabolites.keys():\n",
    "        met2 = met.copy()\n",
    "        if met.id==\"sSUCROSE_b\":\n",
    "            met2.id = \"SUCROSE_ph\"\n",
    "            met = model.metabolites.get_by_id(\"SUCROSE_c\")\n",
    "        elif \"PROTON\" in met.id:\n",
    "            continue\n",
    "        else:\n",
    "            met2.id = met.id.replace(\"_c\",\"_ph\")\n",
    "        met2.compartment = \"ph\"\n",
    "        model.add_metabolites(met2)\n",
    "        \n",
    "        rxn = Reaction(met2.id+\"_exchange\")\n",
    "        rxn.add_metabolites({met2:1})\n",
    "        model.add_reaction(rxn)\n",
    "\n",
    "        rxn = Reaction(met2.id.replace(\"_ph\",\"_phloem_uptake\"),name=met2.id.replace(\"_ph\",\"_phloem_uptake\"))\n",
    "        rxn.add_metabolites({met2:-1,model.metabolites.get_by_id(\"PROTON_e\"):-1,\n",
    "                             met:1,model.metabolites.get_by_id(\"PROTON_c\"):1})\n",
    "        rxn.lower_bound = 0\n",
    "        rxn.upper_bound = 1000\n",
    "        model.add_reaction(rxn)\n",
    "\n",
    "        #print(rxn.reaction)\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "cobra\\core\\metabolite.py:104 UserWarning: 17.612021999304808 is not an integer (in formula C17.612021999304808H32.696359514628455O0.27067590171015216X1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 32.696359514628455 is not an integer (in formula C17.612021999304808H32.696359514628455O0.27067590171015216X1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 0.27067590171015216 is not an integer (in formula C17.612021999304808H32.696359514628455O0.27067590171015216X1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 17.612021999304808 is not an integer (in formula C17.612021999304808H33.180142488275095O2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 33.180142488275095 is not an integer (in formula C17.612021999304808H33.180142488275095O2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 20.612021999304808 is not an integer (in formula C20.612021999304808H38.696359514628455O6.2706759017101525P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 38.696359514628455 is not an integer (in formula C20.612021999304808H38.696359514628455O6.2706759017101525P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 6.2706759017101525 is not an integer (in formula C20.612021999304808H38.696359514628455O6.2706759017101525P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 38.224043998609616 is not an integer (in formula C38.224043998609616H70.39271902925691O6.541351803420305P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 70.39271902925691 is not an integer (in formula C38.224043998609616H70.39271902925691O6.541351803420305P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 6.541351803420305 is not an integer (in formula C38.224043998609616H70.39271902925691O6.541351803420305P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 38.224043998609616 is not an integer (in formula C38.224043998609616H71.39271902925691O3.541351803420305)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 71.39271902925691 is not an integer (in formula C38.224043998609616H71.39271902925691O3.541351803420305)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 3.541351803420305 is not an integer (in formula C38.224043998609616H71.39271902925691O3.541351803420305)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 38.61202199930481 is not an integer (in formula C38.61202199930481H68.1801424882751N7.0O17.000000000000004P3.0S1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 68.1801424882751 is not an integer (in formula C38.61202199930481H68.1801424882751N7.0O17.000000000000004P3.0S1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 17.000000000000004 is not an integer (in formula C38.61202199930481H68.1801424882751N7.0O17.000000000000004P3.0S1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 47.224043998609616 is not an integer (in formula C47.224043998609616H82.39271902925692N3.0O13.541351803420307P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 82.39271902925692 is not an integer (in formula C47.224043998609616H82.39271902925692N3.0O13.541351803420307P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 13.541351803420307 is not an integer (in formula C47.224043998609616H82.39271902925692N3.0O13.541351803420307P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 41.224043998609616 is not an integer (in formula C41.224043998609616H76.39271902925692O11.541351803420305P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 76.39271902925692 is not an integer (in formula C41.224043998609616H76.39271902925692O11.541351803420305P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 11.541351803420305 is not an integer (in formula C41.224043998609616H76.39271902925692O11.541351803420305P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 44.224043998609616 is not an integer (in formula C44.224043998609616H81.39271902925691O8.541351803420305)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 81.39271902925691 is not an integer (in formula C44.224043998609616H81.39271902925691O8.541351803420305)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 8.541351803420305 is not an integer (in formula C44.224043998609616H81.39271902925691O8.541351803420305)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 43.224043998609616 is not an integer (in formula C43.224043998609616H83.39271902925691N1.0O6.541351803420305P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 83.39271902925691 is not an integer (in formula C43.224043998609616H83.39271902925691N1.0O6.541351803420305P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 6.541351803420305 is not an integer (in formula C43.224043998609616H83.39271902925691N1.0O6.541351803420305P1.0)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Biomass flux =39.89663607922134\n"
     ]
    }
   ],
   "source": [
    "from cobra import flux_analysis\n",
    "stem_model = model.copy()\n",
    "\n",
    "k = \"stem\"\n",
    "RXN1 = Reaction(\"Fatty_acid_mix_\"+k)\n",
    "RXN2 = Reaction(\"Fatty_acid_ACP_\"+k)\n",
    "tot = 0\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    #print met\n",
    "    if met in FA:\n",
    "        RXN1.add_metabolites({stem_model.metabolites.get_by_id(met):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "        RXN2.add_metabolites({stem_model.metabolites.get_by_id(FACP[met]):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "        tot = tot+(float(df[df[\"Unnamed: 0\"]==met][k]))\n",
    "if tot==0:\n",
    "    RXN1.add_metabolites({stem_model.metabolites.PALMITATE_p:-1})\n",
    "    RXN2.add_metabolites({stem_model.metabolites.Palmitoyl_ACPs_p:-1})\n",
    "    tot = 1\n",
    "RXN1.add_metabolites({stem_model.metabolites.Fatty_Acids_p:tot})\n",
    "RXN1.lower_bound = 0\n",
    "RXN1.upper_bound = 1000\n",
    "stem_model.add_reaction(RXN1)\n",
    "\n",
    "RXN2.add_metabolites({stem_model.metabolites.Fatty_acyl_ACP_p:tot})\n",
    "RXN2.lower_bound = 0\n",
    "RXN2.upper_bound = 1000\n",
    "stem_model.add_reaction(RXN2)\n",
    "\n",
    "generateMissingFormula(stem_model)\n",
    "\n",
    "rxn = Reaction(\"Biomass_stem_tx\")\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    if met in FA:\n",
    "        continue\n",
    "    else:\n",
    "        rxn.add_metabolites({stem_model.metabolites.get_by_id(met):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "rxn.lower_bound = 0\n",
    "rxn.upper_bound = 1000\n",
    "stem_model.add_reaction(rxn)\n",
    "\n",
    "stem_model = generateStemModel(stem_model)\n",
    "\n",
    "stem_model.reactions.GLC_tx.upper_bound = 0\n",
    "stem_model.reactions.GLC_tx.lower_bound = 0\n",
    "stem_model.reactions.Sucrose_tx.upper_bound = 0\n",
    "stem_model.reactions.Sucrose_tx.lower_bound = 0\n",
    "stem_model.reactions.Photon_tx.upper_bound = 0\n",
    "stem_model.reactions.Photon_tx.lower_bound = 0\n",
    "stem_model.reactions.Biomass_stem_tx.objective_coefficient=1\n",
    "\n",
    "flux_analysis.parsimonious.pfba(stem_model)\n",
    "print(\"Biomass flux =\"+str(stem_model.reactions.Biomass_stem_tx.flux))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'7.659408551696633e-05 CYS_c + 0.6624127272448725 Cellulose_b + 0.0008781839371548684 FRU_c + 0.00010707394063407757 GLC_c + 1.838367529633685e-05 GLY_c + 8.538347701644519e-06 ILE_c + 2.1193398759439072e-05 LEU_c + 9.782907765522723e-06 LYS_c + 5.675545815562533e-06 L_PHOSPHATIDATE_p + 7.975265150758646e-06 MET_c + 2.558895107947629e-05 MYO_INOSITOL_c + 5.93259338961387e-06 PHE_c + 6.861780974383322e-06 PRO_c + 1.1975309152513402e-05 SUC_c + 2.7407527792043203 Starch_b + 7.463115085767901e-05 THR_c + 4.514634314068649e-05 TYR_c + 9.816781301063267e-06 VAL_c + 12.522273984514358 Xylan_b + 4.045921323046277 pALA_b + 0.004964837109152188 pASN_b + 3.1123548220012403 pASP_b + 0.1476511846744479 pCYS_b + 0.17311367516396445 pGLN_b + 2.2107066580017207 pGLU_b + 4.183739547282015 pGLY_b + 2.0708777474801967 pILE_b + 2.893926596282221 pLEU_b + 0.5821794959764395 pLYS_b + 0.5362542739043828 pMET_b + 1.6533728379480637 pPHE_b + 3.042380060845757 pPRO_b + 2.7238288507142943 pSER_b + 1.8980298324803504 pTHR_b + 0.0025742359683561056 pTRP_b + 0.032626784232601035 pTYR_b + 2.7394747885422626 pVAL_b + 0.00013053748895257752 sALA_b + 5.874571538579054e-05 sASP_b + 1.8191461964244266e-05 sCIT_b + 0.0002316400289400987 sFUM_b + 5.9154520918202294e-05 sGABA_b + 8.390288289210637e-05 sGLU_b + 8.154670089457565e-05 sMAL_b + 6.575154972882006e-05 sSER_b + 7.069894496139723e-06 sSUCROSE_b --> '"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stem_model.reactions.Biomass_stem_tx.reaction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generateRootModel(model):\n",
    "    for met in model.reactions.Phloem_output_tx.metabolites.keys():\n",
    "        met2 = met.copy()\n",
    "        if met.id==\"sSUCROSE_b\":\n",
    "            met2.id = \"SUCROSE_ph\"\n",
    "            met = model.metabolites.get_by_id(\"SUCROSE_c\")\n",
    "        elif \"PROTON\" in met.id:\n",
    "            continue\n",
    "        else:\n",
    "            met2.id = met.id.replace(\"_c\",\"_ph\")\n",
    "        met2.compartment = \"ph\"\n",
    "        model.add_metabolites(met2)\n",
    "        \n",
    "        rxn = Reaction(met2.id+\"_exchange\")\n",
    "        rxn.add_metabolites({met2:1})\n",
    "        model.add_reaction(rxn)\n",
    "\n",
    "        rxn = Reaction(met2.id.replace(\"_ph\",\"_phloem_uptake\"),name=met2.id.replace(\"_ph\",\"_phloem_uptake\"))\n",
    "        rxn.add_metabolites({met2:-1,model.metabolites.get_by_id(\"PROTON_e\"):-1,\n",
    "                             met:1,model.metabolites.get_by_id(\"PROTON_c\"):1})\n",
    "        rxn.lower_bound = 0\n",
    "        rxn.upper_bound = 1000\n",
    "        model.add_reaction(rxn)\n",
    "\n",
    "    #add xylem reactions\n",
    "    for met in [\"CAII\",\"MGII\",\"KI\",\"NITRATE\",\"SULFATE\",\"AMMONIUM\",\"WATER\",\"GLT\",\"L_ASPARTATE\",\"ASN\",\"GLN\"]:\n",
    "        met2 = model.metabolites.get_by_id(met+\"_c\").copy()\n",
    "        met2.id = met+\"_xy\"\n",
    "        met2.compartment = \"xy\"\n",
    "\n",
    "        rxn = Reaction(met+\"_exchange\")\n",
    "        rxn.add_metabolites({met2:-1})\n",
    "        rxn.lower_bound = 0\n",
    "        rxn.upper_bound = 1000\n",
    "        model.add_reaction(rxn)\n",
    "\n",
    "        rxn = Reaction(met+\"_xylem_export\")\n",
    "        rxn.add_metabolites({model.metabolites.get_by_id(met+\"_c\"):-1,met2:1})\n",
    "        model.add_reaction(rxn)\n",
    "        \n",
    "    #adding symbiont compartment\n",
    "    from cobra import io\n",
    "    rhizo = io.read_sbml_model(\"Data/iGD1348.xml\")\n",
    "    rhizo.compartments={\"c0\":\"rhizobe cytosol\",\"e0\":\"rhizobe extracellular\"}\n",
    "    model = model+rhizo\n",
    "    \n",
    "    rxn = Reaction(\"Sucrose_exchange_symbiont\")\n",
    "    rxn.name = rxn.id.replace(\"_\",\" \")\n",
    "    rxn.add_metabolites({model.metabolites.get_by_id(\"SUCROSE_c\"):-1,model.metabolites.get_by_id(\"cpd00076[e0]\"):1})\n",
    "    rxn.lower_bound = -1000\n",
    "    rxn.upper_bound = 1000\n",
    "    model.add_reaction(rxn)\n",
    "    \n",
    "    rxn = Reaction(\"Alanine_exchange_symbiont\")\n",
    "    rxn.name = rxn.id.replace(\"_\",\" \")\n",
    "    rxn.add_metabolites({model.metabolites.get_by_id(\"L_ALPHA_ALANINE_c\"):-1,model.metabolites.get_by_id(\"cpd00035[e0]\"):1})\n",
    "    rxn.lower_bound = -1000\n",
    "    rxn.upper_bound = 1000\n",
    "    model.add_reaction(rxn)\n",
    "    \n",
    "    rxn = Reaction(\"Glutamate_exchange_symbiont\")\n",
    "    rxn.name = rxn.id.replace(\"_\",\" \")\n",
    "    rxn.add_metabolites({model.metabolites.get_by_id(\"GLT_c\"):-1,model.metabolites.get_by_id(\"cpd00023[e0]\"):1})\n",
    "    rxn.lower_bound = -1000\n",
    "    rxn.upper_bound = 1000\n",
    "    model.add_reaction(rxn)\n",
    "    \n",
    "    rxn = Reaction(\"Malate_exchange_symbiont\")\n",
    "    rxn.name = rxn.id.replace(\"_\",\" \")\n",
    "    rxn.add_metabolites({model.metabolites.get_by_id(\"MAL_c\"):-1,model.metabolites.get_by_id(\"cpd00130[e0]\"):1})\n",
    "    rxn.lower_bound = -1000\n",
    "    rxn.upper_bound = 1000\n",
    "    model.add_reaction(rxn)\n",
    "    \n",
    "    rxn = Reaction(\"Succinate_exchange_symbiont\")\n",
    "    rxn.name = rxn.id.replace(\"_\",\" \")\n",
    "    rxn.add_metabolites({model.metabolites.get_by_id(\"SUC_c\"):-1,model.metabolites.get_by_id(\"cpd00036[e0]\"):1})\n",
    "    rxn.lower_bound = -1000\n",
    "    rxn.upper_bound = 1000\n",
    "    model.add_reaction(rxn)\n",
    "    \n",
    "    rxn = Reaction(\"Ammonium_exchange_symbiont\")\n",
    "    rxn.name = rxn.id.replace(\"_\",\" \")\n",
    "    rxn.add_metabolites({model.metabolites.get_by_id(\"AMMONIUM_c\"):-1,model.metabolites.get_by_id(\"cpdFixed[c0]\"):1})\n",
    "    rxn.lower_bound = -1000\n",
    "    rxn.upper_bound = 1000\n",
    "    model.add_reaction(rxn)\n",
    "    \n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "cobra\\core\\metabolite.py:104 UserWarning: 17.700746860586182 is not an integer (in formula C17.700746860586182H33.000213348999786O0.3450481094278798X1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 33.000213348999786 is not an integer (in formula C17.700746860586182H33.000213348999786O0.3450481094278798X1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 0.3450481094278798 is not an integer (in formula C17.700746860586182H33.000213348999786O0.3450481094278798X1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 17.700746860586182 is not an integer (in formula C17.700746860586182H33.30951059910088O2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 33.30951059910088 is not an integer (in formula C17.700746860586182H33.30951059910088O2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 20.700746860586182 is not an integer (in formula C20.700746860586182H39.000213348999786O6.3450481094278794P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 39.000213348999786 is not an integer (in formula C20.700746860586182H39.000213348999786O6.3450481094278794P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 6.3450481094278794 is not an integer (in formula C20.700746860586182H39.000213348999786O6.3450481094278794P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 38.401493721172365 is not an integer (in formula C38.401493721172365H71.00042669799957O6.690096218855759P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 71.00042669799957 is not an integer (in formula C38.401493721172365H71.00042669799957O6.690096218855759P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 6.690096218855759 is not an integer (in formula C38.401493721172365H71.00042669799957O6.690096218855759P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 38.401493721172365 is not an integer (in formula C38.401493721172365H72.00042669799957O3.690096218855759)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 72.00042669799957 is not an integer (in formula C38.401493721172365H72.00042669799957O3.690096218855759)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 3.690096218855759 is not an integer (in formula C38.401493721172365H72.00042669799957O3.690096218855759)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 38.70074686058618 is not an integer (in formula C38.70074686058618H68.30951059910089N7.0O17.000000000000004P3.0S1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 68.30951059910089 is not an integer (in formula C38.70074686058618H68.30951059910089N7.0O17.000000000000004P3.0S1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 17.000000000000004 is not an integer (in formula C38.70074686058618H68.30951059910089N7.0O17.000000000000004P3.0S1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 47.401493721172365 is not an integer (in formula C47.401493721172365H83.00042669799959N3.0O13.690096218855759P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 83.00042669799959 is not an integer (in formula C47.401493721172365H83.00042669799959N3.0O13.690096218855759P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 13.690096218855759 is not an integer (in formula C47.401493721172365H83.00042669799959N3.0O13.690096218855759P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 41.401493721172365 is not an integer (in formula C41.401493721172365H77.00042669799959O11.690096218855757P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 77.00042669799959 is not an integer (in formula C41.401493721172365H77.00042669799959O11.690096218855757P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 11.690096218855757 is not an integer (in formula C41.401493721172365H77.00042669799959O11.690096218855757P2.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 44.401493721172365 is not an integer (in formula C44.401493721172365H82.00042669799957O8.690096218855757)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 82.00042669799957 is not an integer (in formula C44.401493721172365H82.00042669799957O8.690096218855757)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 8.690096218855757 is not an integer (in formula C44.401493721172365H82.00042669799957O8.690096218855757)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 43.401493721172365 is not an integer (in formula C43.401493721172365H84.00042669799957N1.0O6.690096218855759P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 84.00042669799957 is not an integer (in formula C43.401493721172365H84.00042669799957N1.0O6.690096218855759P1.0)\n",
      "cobra\\core\\metabolite.py:104 UserWarning: 6.690096218855759 is not an integer (in formula C43.401493721172365H84.00042669799957N1.0O6.690096218855759P1.0)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Biomass flux = 0.0\n"
     ]
    }
   ],
   "source": [
    "from cobra import flux_analysis\n",
    "root_model = model.copy()\n",
    "\n",
    "k = \"root\"\n",
    "RXN1 = Reaction(\"Fatty_acid_mix_\"+k)\n",
    "RXN2 = Reaction(\"Fatty_acid_ACP_\"+k)\n",
    "tot = 0\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    #print met\n",
    "    if met in FA:\n",
    "        RXN1.add_metabolites({root_model.metabolites.get_by_id(met):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "        RXN2.add_metabolites({root_model.metabolites.get_by_id(FACP[met]):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "        tot = tot+(float(df[df[\"Unnamed: 0\"]==met][k]))\n",
    "if tot==0:\n",
    "    RXN1.add_metabolites({root_model.metabolites.PALMITATE_p:-1})\n",
    "    RXN2.add_metabolites({root_model.metabolites.Palmitoyl_ACPs_p:-1})\n",
    "    tot = 1\n",
    "RXN1.add_metabolites({root_model.metabolites.Fatty_Acids_p:tot})\n",
    "RXN1.lower_bound = 0\n",
    "RXN1.upper_bound = 1000\n",
    "root_model.add_reaction(RXN1)\n",
    "    \n",
    "RXN2.add_metabolites({root_model.metabolites.Fatty_acyl_ACP_p:tot})\n",
    "RXN2.lower_bound = 0\n",
    "RXN2.upper_bound = 1000\n",
    "root_model.add_reaction(RXN2)\n",
    "\n",
    "\n",
    "generateMissingFormula(root_model)\n",
    "\n",
    "rxn = Reaction(\"Biomass_root_tx\")\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    if met in FA:\n",
    "        continue\n",
    "    rxn.add_metabolites({root_model.metabolites.get_by_id(met):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "rxn.lower_bound = 0\n",
    "rxn.upper_bound = 1000\n",
    "root_model.add_reaction(rxn)\n",
    "\n",
    "root_model = generateRootModel(root_model)\n",
    "\n",
    "root_model.reactions.GLC_tx.upper_bound = 0\n",
    "root_model.reactions.GLC_tx.lower_bound = 0\n",
    "root_model.reactions.Sucrose_tx.upper_bound = 0\n",
    "root_model.reactions.Sucrose_tx.lower_bound = 0\n",
    "root_model.reactions.Photon_tx.upper_bound = 0\n",
    "root_model.reactions.Photon_tx.lower_bound = 0\n",
    "\n",
    "\n",
    "root_model.reactions.Nitrate_tx.upper_bound = 0\n",
    "root_model.reactions.Nitrate_tx.lower_bound = 0\n",
    "root_model.reactions.NH4_tx.upper_bound = 0\n",
    "root_model.reactions.NH4_tx.lower_bound = 0\n",
    "#root_model.reactions.EX_cpdFixed_e0.upper_bound = 0\n",
    "#root_model.reactions.EX_cpdFixed_e0.lower_bound = 0\n",
    "#root_model.reactions.EX_cpd00013_e0.lower_bound = 0\n",
    "#root_model.reactions.EX_cpd00013_e0.upper_bound = 0 \n",
    "\n",
    "for rxn in root_model.reactions.query(\"EX_cpd\"):\n",
    "    if rxn.id in [#\"EX_cpd00076_e0\",\n",
    "                  \"EX_cpd00001_e0\",\"EX_cpd00067_e0\",\"EX_cpd00011_e0\",\n",
    "                  \"EX_cpd00007_e0\",\"EX_cpd00149_e0\",\"EX_cpd00009_e0\",\n",
    "                  \"EX_cpd00254_e0\",\"EX_cpd00063_e0\",\"EX_cpd00030_e0\",\n",
    "                  \"EX_cpd10515_e0\",\"EX_cpd10516_e0\",\"EX_cpd00034_e0\",\n",
    "                  \"EX_cpd00205_e0\",\"EX_cpd00971_e0\",\"EX_cpd00104_e0\",\n",
    "                  \"EX_cpd00048_e0\",\"EX_cpd11416_e0\",\"EX_cpd00098_e0\",\n",
    "                  \"EX_cpd00528_e0\"]:#\"EX_cpd00182_e0\",\"EX_cpd00311_e0\",\"EX_cpd00064_e0\",\n",
    "        continue\n",
    "    rxn.lower_bound = 0\n",
    "    rxn.upper_bound = 0\n",
    "    \n",
    "for rxn in root_model.reactions.query(\"_phloem_\"):\n",
    "    if \"SUCROSE\" not in rxn.id :\n",
    "        rxn.lower_bound = 0\n",
    "        rxn.upper_bound = 0\n",
    "    else:\n",
    "        rxn.upper_bound = 10\n",
    "        \n",
    "\n",
    "root_model.reactions.EX_cpd00013_e0.lower_bound = 0\n",
    "root_model.reactions.EX_cpd00013_e0.upper_bound = 0\n",
    "root_model.reactions.EX_cpdFixed_e0.upper_bound = 0\n",
    "root_model.reactions.EX_cpdFixed_e0.lower_bound = 0\n",
    "\n",
    "root_model.reactions.GLN_exchange.lower_bound = 0\n",
    "for rxn in root_model.reactions:\n",
    "    rxn.objective_coefficient=0\n",
    "root_model.reactions.Biomass_root_tx.objective_coefficient=1\n",
    "flux_analysis.parsimonious.pfba(root_model)\n",
    "print(\"Biomass flux = \"+str(root_model.reactions.Biomass_root_tx.flux))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "met = root_model.metabolites.get_by_id(\"cpd00013[c0]\")\n",
    "for rxn in met.reactions:\n",
    "    if round(rxn.flux,7)!=0:\n",
    "        eqn = rxn.reaction\n",
    "        for met2 in rxn.metabolites.keys():\n",
    "            eqn = eqn.replace(met2.id,met2.name)\n",
    "        #print(rxn.id+\"\\t\"+str(rxn.metabolites[met]*rxn.flux))\n",
    "        print(rxn.id+\"\\t\"+str(rxn.metabolites[met]*rxn.flux)+\"\\t\"+eqn)\n",
    "        #print(rxn.metabolites)\n",
    "        #print(eqn+\"\\t\"+str(rxn.metabolites[met]))\n",
    "        #print(eqn+\"\\t\"+str(rxn.flux))\n",
    "        #print(\"------------------\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rxn08615\tcpd00387[c0] + cpd15302[c0] <=> cpd00008[c0] + cpd00067[c0] + cpd00155[c0]\t-1000.0\t1000.0\n",
      "ADPglucose_c0 + glycogenn_1_c0 <=> ADP_c0 + H_c0 + Glycogen_c0\n",
      "----------\n",
      "rxn10052\tcpd00018[c0] + cpd00067[c0] + cpd00421[c0] <=> cpd00008[c0] + cpd00012[c0]\t-1000.0\t1000.0\n",
      "AMP_c0 + H_c0 + Triphosphate_c0 <=> ADP_c0 + PPi_c0\n",
      "----------\n",
      "rxn05231\tcpd00008[c0] + cpd00067[c0] + cpd11421[c0] --> cpd00001[c0] + cpd00177[c0] + cpd11420[c0]\t0.0\t1000.0\n",
      "ADP_c0 + H_c0 + trdrd_c0 --> H2O_c0 + dADP_c0 + trdox_c0\n",
      "----------\n"
     ]
    }
   ],
   "source": [
    "for rxn in root_model.metabolites.get_by_id(\"cpd00008[c0]\").reactions:\n",
    "    skip = False\n",
    "    eqn = rxn.reaction\n",
    "    for met2 in rxn.metabolites.keys():\n",
    "        if \"ATP\" in met2.name:\n",
    "            skip = True\n",
    "        eqn = eqn.replace(met2.id,met2.name)\n",
    "    if skip==True:\n",
    "        continue\n",
    "    print(rxn.id+\"\\t\"+rxn.reaction+\"\\t\"+str(rxn.lower_bound)+\"\\t\"+str(rxn.upper_bound))\n",
    "    print(eqn)\n",
    "    print(\"----------\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 H_c0 + N2_c0 --> H2_c0 + 2 Fixed ammonium_c0\n",
      "Biomass flux =0.0\n"
     ]
    }
   ],
   "source": [
    "temp = root_model.copy()\n",
    "#temp.metabolites.get_by_id(\"cpd11620[c0]\").remove_from_model()\n",
    "#temp.metabolites.get_by_id(\"cpd11621[c0]\").remove_from_model()\n",
    "#temp.metabolites.get_by_id(\"cpd00002[c0]\").remove_from_model()\n",
    "\n",
    "#temp.metabolites.get_by_id(\"cpd00008[c0]\").remove_from_model()\n",
    "#temp.metabolites.get_by_id(\"cpd00009[c0]\").remove_from_model()\n",
    "#temp.metabolites.get_by_id(\"cpd00067[e0]\").remove_from_model()\n",
    "#temp.metabolites.get_by_id(\"cpd00067[c0]\").remove_from_model()\n",
    "#temp.metabolites.get_by_id(\"cpd00001[c0]\").remove_from_model()\n",
    "#temp.metabolites.get_by_id(\"cpd00018[c0]\").remove_from_model()\n",
    "\n",
    "rxn = Reaction(\"temp1\")\n",
    "rxn.add_metabolites({temp.metabolites.get_by_id(\"cpdFixed[c0]\"):-1,\n",
    "                     temp.metabolites.get_by_id(\"cpd00013[c0]\"):1})\n",
    "#rxn.add_metabolites({temp.metabolites.get_by_id(\"cpd00013[c0]\"):-1})\n",
    "rxn.upper_bound = 1000\n",
    "rxn.lower_bound = -1000\n",
    "temp.add_reaction(rxn)\n",
    "\n",
    "rxn = Reaction(\"temp2\")\n",
    "rxn.add_metabolites({#temp.metabolites.get_by_id(\"cpd00001[c0]\"):-16,\n",
    "                     #temp.metabolites.get_by_id(\"cpd00002[c0]\"):-16,\n",
    "                     temp.metabolites.get_by_id(\"cpd00528[c0]\"):-1,\n",
    "                     #temp.metabolites.get_by_id(\"cpd11620[c0]\"):-4,\n",
    "                     #temp.metabolites.get_by_id(\"cpd00008[c0]\"):16,\n",
    "                     #temp.metabolites.get_by_id(\"cpd00009[c0]\"):16,\n",
    "                     temp.metabolites.get_by_id(\"cpd00067[c0]\"):-6,\n",
    "                     #temp.metabolites.get_by_id(\"cpd11621[c0]\"):4,\n",
    "                     temp.metabolites.get_by_id(\"cpd11640[c0]\"):1,\n",
    "                     temp.metabolites.get_by_id(\"cpdFixed[c0]\"):2})\n",
    "# # rxn.add_metabolites({temp.metabolites.get_by_id(\"cpd00018[c0]\"):-1})\n",
    "# rxn.lower_bound = -1000\n",
    "temp.add_reaction(rxn)\n",
    "\n",
    "eqn = rxn.reaction\n",
    "for met2 in rxn.metabolites.keys():\n",
    "    eqn = eqn.replace(met2.id,met2.name)\n",
    "print(eqn)\n",
    "\n",
    "flux_analysis.parsimonious.pfba(temp)\n",
    "print(\"Biomass flux =\"+str(temp.reactions.Biomass_root_tx.flux))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "symbiosis_c0\t-0.0\tL_Alanine_c0 + L_Aspartate_c0 + L_Arginine_c0 + 3.0 H_c0 + Glycogen_c0 + Palmitate_c0 + Fixed ammonium_c0 + 0.01 Symbiosis cofactors_c0 --> L_Alanine_e0 + L_Aspartate_e0 + glycogenn_1_c0 + Fixed ammonium_e0\n",
      "cpd00035[c0] + cpd00041[c0] + cpd00051[c0] + 3.0 cpd00067[c0] + cpd00155[c0] + cpd00214[c0] + cpdFixed[c0] + 0.01 cpdSymCoF[c0] --> cpd00035[e0] + cpd00041[e0] + cpd15302[c0] + cpdFixed[e0]\n",
      "temp2\t0.0\t6 H_c0 + N2_c0 --> H2_c0 + 2 Fixed ammonium_c0\n",
      "6 cpd00067[c0] + cpd00528[c0] --> cpd11640[c0] + 2 cpdFixed[c0]\n",
      "rxn06874_c0\t0.0\t16.0 H2O_c0 + 16.0 ATP_c0 + N2_c0 + 4.0 Reducedferredoxin_c0 --> 16.0 ADP_c0 + 16.0 Phosphate_c0 + 6.0 H_c0 + 4.0 Oxidizedferredoxin_c0 + H2_c0 + 2.0 Fixed ammonium_c0\n",
      "16.0 cpd00001[c0] + 16.0 cpd00002[c0] + cpd00528[c0] + 4.0 cpd11620[c0] --> 16.0 cpd00008[c0] + 16.0 cpd00009[c0] + 6.0 cpd00067[c0] + 4.0 cpd11621[c0] + cpd11640[c0] + 2.0 cpdFixed[c0]\n",
      "Ammonium_exchange_symbiont\t0.0\tAMMONIUM <=> Fixed ammonium_c0\n",
      "AMMONIUM_c <=> cpdFixed[c0]\n"
     ]
    }
   ],
   "source": [
    "met = temp.metabolites.get_by_id(\"cpdFixed[c0]\")\n",
    "for rxn in met.reactions:\n",
    "#     if round(rxn.flux,7)!=0:\n",
    "        eqn = rxn.reaction\n",
    "        for met2 in rxn.metabolites.keys():\n",
    "            eqn = eqn.replace(met2.id,met2.name)\n",
    "        #print(rxn.id+\"\\t\"+str(rxn.metabolites[met]*rxn.flux))\n",
    "        print(rxn.id+\"\\t\"+str(rxn.metabolites[met]*rxn.flux)+\"\\t\"+eqn)\n",
    "        print(rxn.reaction)\n",
    "        #print(eqn+\"\\t\"+str(rxn.metabolites[met]))\n",
    "        #print(eqn+\"\\t\"+str(rxn.flux))\n",
    "        #print(\"------------------\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SUCROSE_ph_exchange\t0.7381768150238426\t --> sSUCROSE[b]\n",
      " --> SUCROSE_ph\n",
      "Sucrose_exchange_symbiont\t0.24999999999999398\tSUCROSE <=> Sucrose_e0\n",
      "SUCROSE_c <=> cpd00076[e0]\n",
      "Alanine_exchange_symbiont\t-0.9999999999999749\tL-ALPHA-ALANINE <=> L_Alanine_e0\n",
      "L_ALPHA_ALANINE_c <=> cpd00035[e0]\n"
     ]
    }
   ],
   "source": [
    "for rxn in temp.reactions.query(\"exchange\"):\n",
    "    if round(rxn.flux,7)!=0:\n",
    "        eqn = rxn.reaction\n",
    "        for met2 in rxn.metabolites.keys():\n",
    "            eqn = eqn.replace(met2.id,met2.name)\n",
    "        #print(rxn.id+\"\\t\"+str(rxn.metabolites[met]*rxn.flux))\n",
    "        print(rxn.id+\"\\t\"+str(rxn.flux)+\"\\t\"+eqn)\n",
    "        print(rxn.reaction)\n",
    "        #print(eqn+\"\\t\"+str(rxn.metabolites[met]))\n",
    "        #print(eqn+\"\\t\"+str(rxn.flux))\n",
    "        #print(\"------------------\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H_tx\tPROTON --> \t0.06439159115954657\n",
      "H2O_tx\t <=> WATER\t-2.1479356960515594\n",
      "CO2_tx\t <=> CARBON-DIOXIDE\t-1.5122562605512555\n",
      "O2_tx\t <=> OXYGEN-MOLECULE\t1.2258448279988976\n",
      "Pi_tx\t <=> Pi[e]\t1.3045641055512395e-07\n",
      "SO4_tx\t <=> SULFATE\t0.029312325460245696\n",
      "Biomass_root_tx\t5.8518541628802956e-05 CYS + 0.6480437957395563 Cellulose[b] + 0.0009546732529629342 FRU[c] + 5.117790215895257e-05 GLC + 8.658977494651417e-06 GLY + 1.1282816605744545e-05 ILE + 2.8511982503705806e-05 LEU + 9.375286608625942e-06 LYS + 3.5567809884531068e-06 L-PHOSPHATIDATE + 4.825370511383384e-06 MET + 4.16306145544625e-05 MYO-INOSITOL + 5.569373386168123e-06 PHE + 1.3984135909819174e-05 PRO + 8.356870415782735e-06 SUC + 1.827811330942959 Starch[b] + 5.674989649020359e-05 THR + 2.9858400536810998e-05 TYR + 1.7499479710591047e-05 VAL + 13.015276386519526 Xylan[b] + 3.039393987283876 pALA[b] + 6.898323860987215e-05 pASN[b] + 2.2788993755542424 pASP[b] + 0.11753326040387865 pCYS[b] + 0.22272671636633606 pGLN[b] + 1.7714302375587634 pGLU[b] + 2.989982793945644 pGLY[b] + 1.6819582563230278 pILE[b] + 2.2626803832681324 pLEU[b] + 0.6972499770713101 pLYS[b] + 0.6815785294871021 pMET[b] + 1.4941036190393633 pPHE[b] + 2.2935330780403382 pPRO[b] + 2.1853317275111146 pSER[b] + 1.862190165713056 pTHR[b] + 0.001971534739252036 pTRP[b] + 0.4487363792434114 pTYR[b] + 2.31196934535044 pVAL[b] + 0.00019990306777690502 sALA[b] + 6.389353580619487e-05 sASP[b] + 3.60656309872517e-05 sCIT[b] + 0.00010617565293204374 sFUM[b] + 7.2730968342052e-05 sGABA[b] + 0.00011579419074506035 sGLU[b] + 0.00016809069265177152 sMAL[b] + 7.069956794285573e-05 sSER[b] + 6.222675734205623e-06 sSUCROSE[b] --> \t0.03667822511037485\n",
      "SUCROSE_ph_exchange\t --> sSUCROSE[b]\t0.738176815023844\n",
      "EX_cpd00001_e0\tH2O_e0 <=> \t0.7499999999999826\n",
      "EX_cpd00067_e0\tH_e0 <=> \t0.999999999999975\n",
      "temp1\tNH3_c0 <=> \t-0.9999999999999754\n"
     ]
    }
   ],
   "source": [
    "for rxn in temp.reactions:\n",
    "    if round(rxn.flux,7)==0:\n",
    "        continue\n",
    "    if len(rxn.reactants)==0 or len(rxn.products)==0:\n",
    "        eqn = rxn.reaction\n",
    "        for met2 in rxn.metabolites.keys():\n",
    "            eqn = eqn.replace(met2.id,met2.name)\n",
    "        print(rxn.id+\"\\t\"+eqn+\"\\t\"+str(rxn.flux))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "for met in temp.metabolites:\n",
    "    met.charge = int(met.charge)\n",
    "io.write_sbml_model(temp,\"tempSBML.xml\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['__builtins__',\n",
       " '__cached__',\n",
       " '__doc__',\n",
       " '__file__',\n",
       " '__loader__',\n",
       " '__name__',\n",
       " '__package__',\n",
       " '__path__',\n",
       " '__spec__',\n",
       " 'absolute_import',\n",
       " 'dict',\n",
       " 'from_json',\n",
       " 'from_yaml',\n",
       " 'json',\n",
       " 'load_json_model',\n",
       " 'load_matlab_model',\n",
       " 'load_yaml_model',\n",
       " 'mat',\n",
       " 'model_from_dict',\n",
       " 'model_to_dict',\n",
       " 'read_sbml_model',\n",
       " 'save_json_model',\n",
       " 'save_matlab_model',\n",
       " 'save_yaml_model',\n",
       " 'sbml',\n",
       " 'to_json',\n",
       " 'to_yaml',\n",
       " 'validate_sbml_model',\n",
       " 'write_sbml_model',\n",
       " 'yaml']"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dir(io)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra import flux_analysis\n",
    "seed_model = model.copy()\n",
    "\n",
    "k = \"seed\"\n",
    "RXN1 = Reaction(\"Fatty_acid_mix_\"+k)\n",
    "RXN2 = Reaction(\"Fatty_acid_ACP_\"+k)\n",
    "tot = 0\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    #print met\n",
    "    if met in FA:\n",
    "        RXN1.add_metabolites({seed_model.metabolites.get_by_id(met):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "        RXN2.add_metabolites({seed_model.metabolites.get_by_id(FACP[met]):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "        tot = tot+(float(df[df[\"Unnamed: 0\"]==met][k]))\n",
    "if tot==0:\n",
    "    RXN1.add_metabolites({seed_model.metabolites.PALMITATE_p:-1})\n",
    "    RXN2.add_metabolites({seed_model.metabolites.Palmitoyl_ACPs_p:-1})\n",
    "    tot = 1\n",
    "RXN1.add_metabolites({seed_model.metabolites.Fatty_Acids_p:tot})\n",
    "RXN1.lower_bound = 0\n",
    "RXN1.upper_bound = 1000\n",
    "seed_model.add_reaction(RXN1)\n",
    "    \n",
    "RXN2.add_metabolites({seed_model.metabolites.Fatty_acyl_ACP_p:tot})\n",
    "RXN2.lower_bound = 0\n",
    "RXN2.upper_bound = 1000\n",
    "seed_model.add_reaction(RXN2)\n",
    "\n",
    "\n",
    "generateMissingFormula(seed_model)\n",
    "\n",
    "rxn = Reaction(\"Biomass_seed_tx\")\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    if met in FA:\n",
    "        continue\n",
    "    rxn.add_metabolites({seed_model.metabolites.get_by_id(met):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "rxn.lower_bound = 0\n",
    "rxn.upper_bound = 1000\n",
    "seed_model.add_reaction(rxn)\n",
    "\n",
    "seed_model.reactions.GLC_tx.upper_bound = 0\n",
    "seed_model.reactions.GLC_tx.lower_bound = 0\n",
    "seed_model.reactions.Sucrose_tx.upper_bound = 1\n",
    "seed_model.reactions.Sucrose_tx.lower_bound = 1\n",
    "seed_model.reactions.Photon_tx.upper_bound = 0\n",
    "seed_model.reactions.Photon_tx.lower_bound = 0\n",
    "seed_model.reactions.Biomass_seed_tx.objective_coefficient=1\n",
    "sol=flux_analysis.parsimonious.pfba(seed_model)\n",
    "print(\"Biomass flux =\"+str(seed_model.reactions.Biomass_seed_tx.flux))\n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_model.reactions.Biomass_seed_tx.reaction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate diel leaf model\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra import flux_analysis\n",
    "leaf_model = model.copy()\n",
    "\n",
    "k = \"leaf\"\n",
    "RXN1 = Reaction(\"Fatty_acid_mix_\"+k)\n",
    "RXN2 = Reaction(\"Fatty_acid_ACP_\"+k)\n",
    "tot = 0\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    #print met\n",
    "    if met in FA:\n",
    "        RXN1.add_metabolites({leaf_model.metabolites.get_by_id(met):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "        RXN2.add_metabolites({leaf_model.metabolites.get_by_id(FACP[met]):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "        tot = tot+(float(df[df[\"Unnamed: 0\"]==met][k]))\n",
    "if tot==0:\n",
    "    RXN1.add_metabolites({leaf_model.metabolites.PALMITATE_p:-1})\n",
    "    RXN2.add_metabolites({leaf_model.metabolites.Palmitoyl_ACPs_p:-1})\n",
    "    tot = 1\n",
    "RXN1.add_metabolites({leaf_model.metabolites.Fatty_Acids_p:tot})\n",
    "RXN1.lower_bound = 0\n",
    "RXN1.upper_bound = 1000\n",
    "leaf_model.add_reaction(RXN1)\n",
    "    \n",
    "RXN2.add_metabolites({leaf_model.metabolites.Fatty_acyl_ACP_p:tot})\n",
    "RXN2.lower_bound = 0\n",
    "RXN2.upper_bound = 1000\n",
    "leaf_model.add_reaction(RXN2)\n",
    "\n",
    "\n",
    "generateMissingFormula(leaf_model)\n",
    "\n",
    "rxn = Reaction(\"Biomass_leaf_tx\")\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    if met in FA:\n",
    "        continue\n",
    "    rxn.add_metabolites({leaf_model.metabolites.get_by_id(met):-1*float(df[df[\"Unnamed: 0\"]==met][k])})\n",
    "\n",
    "rxn.lower_bound = 0\n",
    "rxn.upper_bound = 1000\n",
    "leaf_model.add_reaction(rxn)\n",
    "\n",
    "leaf_model.reactions.GLC_tx.upper_bound = 0\n",
    "leaf_model.reactions.GLC_tx.lower_bound = 0\n",
    "leaf_model.reactions.Sucrose_tx.upper_bound = 1\n",
    "leaf_model.reactions.Sucrose_tx.lower_bound = 1\n",
    "leaf_model.reactions.Photon_tx.upper_bound = 0\n",
    "leaf_model.reactions.Photon_tx.lower_bound = 0\n",
    "leaf_model.reactions.Biomass_leaf_tx.objective_coefficient=1\n",
    "sol=flux_analysis.parsimonious.pfba(leaf_model)\n",
    "print(\"Biomass flux =\"+str(leaf_model.reactions.Biomass_leaf_tx.flux))\n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from sweetlovegroup.transform import setupC3DielModel\n",
    "leaf_model = setupC3DielModel(leaf_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "leaf_model.reactions.Phloem_output_tx1.objective_coefficient = 0\n",
    "leaf_model.reactions.Phloem_output_tx2.objective_coefficient = 0\n",
    "\n",
    "leaf_model.reactions.Photon_tx1.upper_bound = 500\n",
    "leaf_model.reactions.Photon_tx1.lower_bound = 500\n",
    "\n",
    "leaf_model.reactions.Biomass_leaf_tx1.objective_coefficient = 0\n",
    "\n",
    "from cobra.flux_analysis import pfba\n",
    "sol = pfba(leaf_model)\n",
    "print(\"Phloem export rate =\"+str(sol.fluxes[\"diel_biomass\"]*4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "growing_leaf_model = leaf_model.copy()\n",
    "met1=Metabolite(\"Leaf_biomass1\")\n",
    "rxn = growing_leaf_model.reactions.Biomass_leaf_tx1\n",
    "rxn.add_metabolites({met1:1})\n",
    "\n",
    "met2=Metabolite(\"Leaf_biomass2\")\n",
    "rxn = growing_leaf_model.reactions.Biomass_leaf_tx2\n",
    "rxn.add_metabolites({met2:1})\n",
    "\n",
    "growing_leaf_model.reactions.diel_biomass.add_metabolites({growing_leaf_model.metabolites.X_Phloem_contribution_t1:3,\n",
    "                                                           growing_leaf_model.metabolites.X_Phloem_contribution_t2:1,\n",
    "                                                           met1:-3, met2:-1})\n",
    "\n",
    "from cobra.flux_analysis import pfba\n",
    "sol = pfba(growing_leaf_model)\n",
    "print(\"Phloem export rate =\"+str(sol.fluxes[\"diel_biomass\"]*4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}