Generate MAIZE tissue specific models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modelling metabolism in MAIZE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating MAIZE 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 = \"./../plantcoremetabolism-model/PlantCoreMetabolism_v2_1_0.xml\"\n",
    "model = io.read_sbml_model(fname)\n",
    "model.solver='glpk'\n",
    "display.clear_output()\n",
    "\n",
    "\n",
    "\n",
    "#change Biomass_tx to Protein_biomass\n",
    "rxn = model.reactions.get_by_id(\"AraCore_Biomass_tx\")\n",
    "rxn.remove_from_model()\n",
    "\n",
    "#remove charged state of PGPs since they are minor (0.03)\n",
    "model=removeSpecificMetChargedState(model,[\"aL_1_PHOSPHATIDYL_GLYCEROL_P_p\",\"aL_1_PHOSPHATIDYL_GLYCEROL_P_m\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_csv(\"Data/biomass_maize.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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "biomass = df.set_index(\"Unnamed: 0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use biomass composition to generate biomass equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run stem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "Biomass flux =0.8937627272359395\n"
     ]
    }
   ],
   "source": [
    "from cobra import flux_analysis\n",
    "stem_model = model.copy()\n",
    "\n",
    "#change Biomass_tx to Protein_biomass\n",
    "rxn = stem_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  = stem_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",
    "k = \"stem\"\n",
    "\n",
    "stem_model = updateFAcomposition(stem_model,k,biomass)\n",
    "\n",
    "rxn = Reaction(\"Biomass_stem_tx\")\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    if met in FA or float(df[df[\"Unnamed: 0\"]==met][k])==0:\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",
    "for rxn in stem_model.reactions:\n",
    "    if rxn.id==\"SUCROSE_ph_exchange\":\n",
    "        rxn.lower_bound = 1\n",
    "        rxn.upper_bound = 1\n",
    "    elif \"ph_exchange\" in rxn.id:\n",
    "        rxn.lower_bound = 0\n",
    "        rxn.upper_bound = 0\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": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.6686247303468177 Cellulose_b + 0.5082440106752529 GLC_c + 0.08059221386706178 Protein_b + 0.02923548907713272 Starch_b + 0.2814303970207669 Xylan_b + 0.21817790431973 sSUCROSE_b --> '"
      ]
     },
     "execution_count": 6,
     "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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "Biomass flux = 230.2209885659986\n"
     ]
    }
   ],
   "source": [
    "from cobra import flux_analysis\n",
    "root_model = model.copy()\n",
    "\n",
    "#change Biomass_tx to Protein_biomass\n",
    "rxn = root_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  = root_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",
    "k = \"root\"\n",
    "root_model = updateFAcomposition(root_model,k,biomass)\n",
    "\n",
    "rxn = Reaction(\"Biomass_root_tx\")\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    if met in FA or float(df[df[\"Unnamed: 0\"]==met][k])==0:\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.NH4_tx.upper_bound = 0\n",
    "root_model.reactions.NH4_tx.lower_bound = 0\n",
    " \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",
    "\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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.6686247303468177 Cellulose_b + 0.5082440106752529 GLC_c + 0.08059221386706178 Protein_b + 0.02923548907713272 Starch_b + 0.2814303970207669 Xylan_b + 0.21817790431973 sSUCROSE_b --> '"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stem_model.reactions.Biomass_stem_tx.reaction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.07873630218930545\n",
      "Biomass flux =0.27622051577148315\n"
     ]
    }
   ],
   "source": [
    "from cobra import flux_analysis\n",
    "seed_model = model.copy()\n",
    "\n",
    "#change Biomass_tx to Protein_biomass\n",
    "rxn = seed_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  = seed_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",
    "seed_model = generateSeedModel(seed_model)\n",
    "\n",
    "k = \"seed\"\n",
    "seed_model = updateFAcomposition(seed_model,k,biomass)\n",
    "\n",
    "\n",
    "for rxn in seed_model.reactions:\n",
    "    if rxn.id==\"SUCROSE_ph_exchange\":\n",
    "        rxn.lower_bound = 1\n",
    "        rxn.upper_bound = 1\n",
    "    elif \"ph_exchange\" in rxn.id:\n",
    "        rxn.lower_bound = 0\n",
    "        rxn.upper_bound = 0\n",
    "\n",
    "\n",
    "\n",
    "rxn = Reaction(\"Biomass_seed_tx\")\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    if met in FA or float(df[df[\"Unnamed: 0\"]==met][k])==0:\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 = 0\n",
    "seed_model.reactions.Sucrose_tx.lower_bound = 0\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": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.02611208360102965 ARG_c + 0.008934304834262224 COUMARATE_c + 0.028509562182814783 CYS_c + 0.7288080463514958 Cellulose_b + 0.04677142993318231 GLY_c + 0.017496442114094343 HIS_c + 0.03545891941720899 ILE_c + 0.1290514502197219 LEU_c + 0.019058176219910104 LYS_c + 0.010598781120344694 MET_c + 0.038073359759337434 PHE_c + 0.09372941887934014 PRO_c + 0.2351385482160249 Protein_b + 4.512261602090815 Starch_b + 0.03321020919071476 THR_c + 0.002891654256889925 TRP_c + 0.02558300443406753 TYR_c + 0.0453178430025594 VAL_c + 0.13127251840276788 sALA_b + 0.05786093375103354 sASP_b + 0.1307736574268864 sGLU_b + 0.059184588589518985 sSER_b --> '"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seed_model.reactions.Biomass_seed_tx.reaction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate diel leaf model\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def generateMultiphaseModel(model,tags):\n",
    "    models=dict()\n",
    "    for i in range(0,len(tags)):\n",
    "        temp=model.copy()\n",
    "        for met in temp.metabolites:\n",
    "            met.id = met.id+str(i+1)\n",
    "            met.compartment = met.compartment+str(i+1)\n",
    "            met.name = met.name+\" \"+tags[i]\n",
    "        for rxn in temp.reactions:\n",
    "            rxn.id = rxn.id+str(i+1)\n",
    "            rxn.name = rxn.name+\" \"+tags[i]\n",
    "        if i==0:\n",
    "            modelOut = temp.copy()\n",
    "        else:\n",
    "            modelOut = modelOut+temp\n",
    "            for met in temp.metabolites:\n",
    "                if met.id not in modelOut.metabolites:\n",
    "                    met2 = met.copy()\n",
    "                    #print(met2.id)\n",
    "                    modelOut.add_metabolites(met2)\n",
    "        for comp in temp.compartments:\n",
    "            modelOut.compartments={str(comp):model.compartments[comp[0:len(comp)-1]]+\" \"+tags[i]}\n",
    "    return modelOut\n",
    "\n",
    "\n",
    "def addLinkers(model,baseMets,from_phase,to_phase,reversible=True):\n",
    "    from cobra.core import Reaction\n",
    "    for met in baseMets.keys():\n",
    "        baseMet = met\n",
    "    rxn = Reaction(baseMet+\"_linker\"+from_phase+to_phase)\n",
    "    rxn.name = rxn.id\n",
    "    baseMets = baseMets[baseMet]\n",
    "    for baseMetID in baseMets.keys():\n",
    "        met1 = model.metabolites.get_by_id(baseMetID+from_phase)\n",
    "        met2 = model.metabolites.get_by_id(baseMetID+to_phase)\n",
    "        rxn.add_metabolites({met1:-1*baseMets[baseMetID],met2:1*baseMets[baseMetID]})\n",
    "    rxn.upper_bound = 1000\n",
    "    if reversible:\n",
    "        rxn.lower_bound = -1000\n",
    "    else:\n",
    "        rxn.lower_bound = 0\n",
    "    model.add_reaction(rxn)\n",
    "    return model\n",
    "    \n",
    "\n",
    "def identifyFractionalMets(model):\n",
    "    import re\n",
    "    fractionMets=dict()\n",
    "    for rxn in model.reactions:\n",
    "        for met in rxn.metabolites.keys():\n",
    "            prefix=\"\"\n",
    "            a=re.search(\"^a{1,3}\",met.id)\n",
    "            anion=\"\"\n",
    "            if a:\n",
    "                anion=a.group(0)\n",
    "                prefix=anion\n",
    "            b=re.search(\"^b{1,3}\",met.id)\n",
    "            basic=\"\"\n",
    "            if b:\n",
    "                basic=b.group(0)\n",
    "                prefix=basic\n",
    "            if (not prefix == \"\"):\n",
    "                mainMet=model.metabolites.get_by_id(met.id[len(prefix):])\n",
    "                fractionMets[met.id]=mainMet.id\n",
    "                #fractionMets[met.id]=prefix\n",
    "    return fractionMets\n",
    "\n",
    "\n",
    "# def mergeCompartments(model,comp1,comp2):\n",
    "#     temp = model.copy()\n",
    "#     for met in model.metabolites:\n",
    "#         if met.compartment==comp2:\n",
    "#             M_orig = met.id.replace(comp2,comp1)\n",
    "#             if M_orig in model.metabolites:\n",
    "#                 M_orig = model.metabolites.get_by_id(M_orig)\n",
    "#                 for rxn in M_orig.reactions:\n",
    "#                     M2 = model.metabolites.get_by_id(M_orig.id)\n",
    "#                     coeff = rxn.metabolites[M_orig]\n",
    "#                     rxn.add_metabolites({M_orig:-1*coeff,M2:coeff})\n",
    "                    \n",
    "def addUpperAndLowerBoundsFromFile(model,filename):\n",
    "    import pandas as pd\n",
    "    df = pd.read_csv(filename,sep=\"\\t\")\n",
    "    \n",
    "    for rxn in df[\"reaction\"]:\n",
    "        rxn = model.reactions.get_by_id(str(rxn))\n",
    "        rxn.lower_bound = float(df[df[\"reaction\"]==rxn.id][\"lower-bound\"])\n",
    "        rxn.upper_bound = float(df[df[\"reaction\"]==rxn.id][\"upper-bound\"])\n",
    "    \n",
    "    return model\n",
    "\n",
    "\n",
    "def generateC4leafModel(model,constraintsFile):\n",
    "    tags = [\"mesophyll day\",\"mesophyll night\",\"bundlesheath day\",\"bundlesheath night\"]\n",
    "    C4_model = generateMultiphaseModel(model,tags)\n",
    "\n",
    "    mets1 = [\"STARCH_p\",\"MAL_v\",\"FUM_v\",\"CIT_v\",\"NITRATE_v\",\"SUCROSE_v\",\"GLC_v\",\"FRU_v\"]\n",
    "    for met in mets1:\n",
    "        if met == \"MAL_v\":\n",
    "            mets = {met:{\"MAL_v\":0.7,\"aMAL_v\":0.3}}\n",
    "        elif met == \"FUM_v\":\n",
    "            mets = {met:{\"FUM_v\":0.92,\"aFUM_v\":0.08}}\n",
    "        elif met == \"CIT_v\":\n",
    "            mets = {met:{\"CIT_v\":0.5,\"aCIT_v\":0.5}}\n",
    "        else:\n",
    "            mets = {met:{met:1}}\n",
    "        C4_model = addLinkers(C4_model,mets,\"1\",\"2\")\n",
    "        C4_model = addLinkers(C4_model,mets,\"3\",\"4\")\n",
    "\n",
    "    mets2 = [\"ASN_v\",\"GLT_v\",\"L_ALPHA_ALANINE_v\",\"LEU_v\",\"bHIS_v\",\"MET_v\",\"LYS_v\",\"THR_v\",\"VAL_v\",\"PHE_v\",\n",
    "            \"4_AMINO_BUTYRATE_v\",\"SER_v\",\"ARG_v\",\"GLN_v\",\"GLY_v\",\"L_ASPARTATE_v\",\"ILE_v\",\"PRO_v\",\"CYS_v\",\n",
    "            \"TRP_v\",\"TYR_v\"]\n",
    "    for met in mets2:\n",
    "        mets = {met:{met:1}}\n",
    "        C4_model = addLinkers(C4_model,mets,\"1\",\"2\",reversible=False)\n",
    "        C4_model = addLinkers(C4_model,mets,\"3\",\"4\",reversible=False)\n",
    "\n",
    "\n",
    "    mets3 = [\"MAL_c\",\"PYRUVATE_c\",\"ASN_c\",\"GLT_c\",\"L_ALPHA_ALANINE_c\",\"LEU_c\",\n",
    "            \"MET_c\",\"LYS_c\",\"HIS_c\",\"THR_c\",\"VAL_c\",\"PHE_c\",\"4_AMINO_BUTYRATE_c\",\"SER_c\",\"ARG_c\",\"GLN_c\",\n",
    "            \"GLY_c\",\"L_ASPARTATE_c\",\"ILE_c\",\"PRO_c\",\"CYS_c\",\"TRP_c\",\"TYR_c\",\"G3P_c\",\"GAP_c\",\"Pi_c\",\n",
    "            \"PHOSPHO_ENOL_PYRUVATE_c\"] #\"SUCROSE_c\",\"MALTOSE_p\",\n",
    "    for met in mets3:\n",
    "        if met==\"Pi_c\":\n",
    "            mets = {met:{\"Pi_c\":0.7,\"aPi_c\":0.3}}\n",
    "        else:\n",
    "            mets = {met:{met:1}}\n",
    "        C4_model = addLinkers(C4_model,mets,\"1\",\"3\")\n",
    "        C4_model = addLinkers(C4_model,mets,\"2\",\"4\")\n",
    "\n",
    "\n",
    "    C4_model = addUpperAndLowerBoundsFromFile(C4_model,constraintsFile)   \n",
    "\n",
    "    #Nitrate ratio day:night\n",
    "    met = Metabolite(\"Nitrate_pseudo\")\n",
    "    met.name = \"Nitrate_pseudo\"\n",
    "    met.compartment = \"f\"\n",
    "    met.charge = 0\n",
    "    for i in (1,3):\n",
    "        C4_model.reactions.get_by_id(str(\"Nitrate_tx\" + str(i))).add_metabolites({met : 2})\n",
    "    for i in (2,4):\n",
    "        C4_model.reactions.get_by_id(str(\"Nitrate_tx\" + str(i))).add_metabolites({met : -3})\n",
    "\n",
    "    #Rubisco oxygenase\n",
    "    for i in range(1,5):\n",
    "        met = Metabolite(\"Rubisco_pseudo\" + str(i))\n",
    "        met.name = \"Rubisco_pseudo\" + str(i)\n",
    "        met.compartment = \"f\"\n",
    "        met.charge = 0\n",
    "        C4_model.reactions.get_by_id(str(\"RIBULOSE_BISPHOSPHATE_CARBOXYLASE_RXN_p\" + str(i))).add_metabolites({met : 1})\n",
    "        C4_model.reactions.get_by_id(str(\"RXN_961_p\" + str(i))).add_metabolites({met : -20})\n",
    "\n",
    "\n",
    "    #Maintainence costs\n",
    "    for i in range(1,5):\n",
    "        C4_model.reactions.get_by_id(str(\"ATPase_tx\" + str(i))).upper_bound = C4_model.reactions.get_by_id(str(\"ATPase_tx\" + str(i))).lower_bound = 19.5\n",
    "        met = Metabolite(str(\"Maint_pseudo\" + str(i)))\n",
    "        met.name = str(\"Maint_pseudo\" + str(i))\n",
    "        met.compartment = \"f\"\n",
    "        met.charge = 0\n",
    "        C4_model.reactions.get_by_id(str(\"ATPase_tx\" + str(i))).add_metabolites({met : 1})\n",
    "        C4_model.reactions.get_by_id(str(\"NADPHoxc_tx\" + str(i))).add_metabolites({met : -3})\n",
    "        C4_model.reactions.get_by_id(str(\"NADPHoxm_tx\" + str(i))).add_metabolites({met : -3})\n",
    "        C4_model.reactions.get_by_id(str(\"NADPHoxp_tx\" + str(i))).add_metabolites({met : -3})\n",
    "\n",
    "\n",
    "\n",
    "    for i in range(3,5):\n",
    "        for met in model.metabolites:\n",
    "            if \"CARBON_DIOXIDE\" in met.id or \"OXYGEN_MOLECULE\" in met.id or \"Photon\" in met.id:\n",
    "                continue\n",
    "            if met.compartment == \"e\"+str(i):\n",
    "                rxn = Reaction(met.id+\"_e\"+str(i-2)+str(i))\n",
    "                rxn.add_metabolites({met:1,C4_model.metabolites.get_by_id(met.id.replace(met.compartment,\"e\"+str(i-2))):-1})\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    met1 = Metabolite(\"X_Phloem_contribution_day\",name=\"Phloem output during the day\",compartment=\"b1\")\n",
    "    C4_model.reactions.get_by_id(\"Phloem_output_tx1\").remove_from_model()\n",
    "    C4_model.reactions.get_by_id(\"Phloem_output_tx3\").add_metabolites({met1:1})\n",
    "    met2 = Metabolite(\"X_Phloem_contribution_night\",name=\"Phloem output during at night\",compartment=\"b1\")\n",
    "    C4_model.reactions.get_by_id(\"Phloem_output_tx2\").remove_from_model()\n",
    "    C4_model.reactions.get_by_id(\"Phloem_output_tx4\").add_metabolites({met2:1})\n",
    "\n",
    "    rxn = Reaction(\"diel_biomass\")\n",
    "    rxn.add_metabolites({met1:-0.75,met2:-0.25})\n",
    "    rxn.lower_bound = 0\n",
    "    rxn.upper_bound = 1000\n",
    "    C4_model.add_reaction(rxn)\n",
    "    rxn.objective_coefficient=1\n",
    "    return C4_model\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "Biomass flux =0.0\n"
     ]
    }
   ],
   "source": [
    "from cobra import flux_analysis\n",
    "leaf_model = model.copy()\n",
    "\n",
    "#change Biomass_tx to Protein_biomass\n",
    "rxn = leaf_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  = leaf_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",
    "k = \"leaf\"\n",
    "leaf_model = updateFAcomposition(leaf_model,k,biomass)\n",
    "\n",
    "rxn = Reaction(\"Biomass_leaf_tx\")\n",
    "for met in df[\"Unnamed: 0\"]:\n",
    "    if met in FA or float(df[df[\"Unnamed: 0\"]==met][k])==0:\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 = 0\n",
    "leaf_model.reactions.Sucrose_tx.lower_bound = 0\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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Phloem output = \n",
      "5.054115122057934\n"
     ]
    }
   ],
   "source": [
    "leaf_model_source = generateC4leafModel(leaf_model,\"Maize_Leaf_Constraints.txt\")\n",
    "\n",
    "leaf_model_source.reactions.Photon_tx1.upper_bound = 600\n",
    "leaf_model_source.reactions.Photon_tx3.upper_bound = 400\n",
    "\n",
    "sol = flux_analysis.parsimonious.pfba(leaf_model_source)\n",
    "print(\"Phloem output = \")\n",
    "print(leaf_model_source.reactions.diel_biomass.flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CO2_tx1 flux = 51.187972893722495\n",
      "O2_tx1 flux = -71.75000000000001\n",
      "PSII_RXN_p1 flux = 75.00000000000001\n",
      "RIBULOSE_BISPHOSPHATE_CARBOXYLASE_RXN_p1 flux = 0.0\n",
      "PEPCARBOXYKIN_RXN_c1 flux = 0.0\n",
      "1_PERIOD_1_PERIOD_1_PERIOD_39_RXN_m1 flux = 0.0\n",
      "MALIC_NADP_RXN_p1 flux = 0.0\n",
      "MALIC_NADP_RXN_c1 flux = 0.0\n",
      "------\n",
      "CO2_tx2 flux = -0.2541591153634801\n",
      "O2_tx2 flux = 4.3177860449007\n",
      "PSII_RXN_p2 flux = 0.0\n",
      "RIBULOSE_BISPHOSPHATE_CARBOXYLASE_RXN_p2 flux = 0.0\n",
      "PEPCARBOXYKIN_RXN_c2 flux = 0.0\n",
      "1_PERIOD_1_PERIOD_1_PERIOD_39_RXN_m2 flux = 0.0\n",
      "MALIC_NADP_RXN_p2 flux = 0.0\n",
      "MALIC_NADP_RXN_c2 flux = 0.0\n",
      "------\n",
      "CO2_tx3 flux = 0.0\n",
      "O2_tx3 flux = 8.05975887139129\n",
      "PSII_RXN_p3 flux = 0.0\n",
      "RIBULOSE_BISPHOSPHATE_CARBOXYLASE_RXN_p3 flux = 64.13011828521721\n",
      "PEPCARBOXYKIN_RXN_c3 flux = 0.0\n",
      "1_PERIOD_1_PERIOD_1_PERIOD_39_RXN_m3 flux = 0.0\n",
      "MALIC_NADP_RXN_p3 flux = 62.286155054238414\n",
      "MALIC_NADP_RXN_c3 flux = 0.0\n",
      "------\n",
      "CO2_tx4 flux = 0.0\n",
      "O2_tx4 flux = 7.376547970812598\n",
      "PSII_RXN_p4 flux = 0.0\n",
      "RIBULOSE_BISPHOSPHATE_CARBOXYLASE_RXN_p4 flux = 0.0\n",
      "PEPCARBOXYKIN_RXN_c4 flux = 0.0\n",
      "1_PERIOD_1_PERIOD_1_PERIOD_39_RXN_m4 flux = 0.0\n",
      "MALIC_NADP_RXN_p4 flux = 0.0\n",
      "MALIC_NADP_RXN_c4 flux = 0.0\n",
      "------\n"
     ]
    }
   ],
   "source": [
    "for i in range(1,5):\n",
    "    for rxn in [\"CO2_tx\",\"O2_tx\",\"PSII_RXN_p\",\"RIBULOSE_BISPHOSPHATE_CARBOXYLASE_RXN_p\",\n",
    "                \"PEPCARBOXYKIN_RXN_c\",\"1_PERIOD_1_PERIOD_1_PERIOD_39_RXN_m\",\"MALIC_NADP_RXN_p\",\n",
    "                \"MALIC_NADP_RXN_c\"]:\n",
    "        print(rxn+str(i)+\" flux = \"+str(sol.fluxes[rxn+str(i)]))\n",
    "    print(\"------\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Biomass accumulation rate rate =0.0\n"
     ]
    }
   ],
   "source": [
    "growing_leaf_model = leaf_model_source.copy()\n",
    "\n",
    "growing_leaf_model.reactions.Biomass_leaf_tx1.upper_bound=1000\n",
    "growing_leaf_model.reactions.Biomass_leaf_tx2.upper_bound=1000\n",
    "growing_leaf_model.reactions.Biomass_leaf_tx3.upper_bound=1000\n",
    "growing_leaf_model.reactions.Biomass_leaf_tx4.upper_bound=1000\n",
    "\n",
    "met1=Metabolite(\"Leaf_biomass1\",compartment=\"c1\")\n",
    "rxn = growing_leaf_model.reactions.Biomass_leaf_tx1\n",
    "rxn.add_metabolites({met1:1})\n",
    "\n",
    "met2=Metabolite(\"Leaf_biomass2\",compartment=\"c2\")\n",
    "rxn = growing_leaf_model.reactions.Biomass_leaf_tx2\n",
    "rxn.add_metabolites({met2:1})\n",
    "\n",
    "met3=Metabolite(\"Leaf_biomass3\",compartment=\"c3\")\n",
    "rxn = growing_leaf_model.reactions.Biomass_leaf_tx3\n",
    "rxn.add_metabolites({met3:1})\n",
    "\n",
    "met4=Metabolite(\"Leaf_biomass4\",compartment=\"c4\")\n",
    "rxn = growing_leaf_model.reactions.Biomass_leaf_tx4\n",
    "rxn.add_metabolites({met4:1})\n",
    "\n",
    "growing_leaf_model.reactions.diel_biomass.add_metabolites({growing_leaf_model.metabolites.X_Phloem_contribution_day:3,\n",
    "                                                           growing_leaf_model.metabolites.X_Phloem_contribution_night:1,\n",
    "                                                           met1:-0.75, met2:-0.25,met3:-0.75, met4:-0.25})\n",
    "\n",
    "\n",
    "from cobra.flux_analysis import pfba\n",
    "sol = pfba(growing_leaf_model)\n",
    "print(\"Biomass accumulation rate rate =\"+str(sol.fluxes[\"diel_biomass\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Protein_b\n",
      "-0.09299016640000002\n",
      "Protein_b\n",
      "-0.09299016640000002\n",
      "Protein_b1\n",
      "-0.09299016640000002\n",
      "Protein_b2\n",
      "-0.09299016640000002\n",
      "Protein_b3\n",
      "-0.09299016640000002\n",
      "Protein_b4\n",
      "-0.09299016640000002\n",
      "X_Phloem_contribution_day\n",
      "None\n",
      "X_Phloem_contribution_night\n",
      "None\n",
      "Protein_b1\n",
      "-0.09299016640000002\n",
      "Protein_b2\n",
      "-0.09299016640000002\n",
      "Protein_b3\n",
      "-0.09299016640000002\n",
      "Protein_b4\n",
      "-0.09299016640000002\n",
      "X_Phloem_contribution_day\n",
      "None\n",
      "X_Phloem_contribution_night\n",
      "None\n",
      "Leaf_biomass1\n",
      "None\n",
      "Leaf_biomass2\n",
      "None\n",
      "Leaf_biomass3\n",
      "None\n",
      "Leaf_biomass4\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "for Tmodel in [root_model,stem_model,leaf_model_source,growing_leaf_model]:\n",
    "    for met in Tmodel.metabolites:\n",
    "        if type(met.charge)!=int:\n",
    "            print(met.id)\n",
    "            print(met.charge)\n",
    "            met.charge=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "io.write_sbml_model(root_model,\"C:\\\\Users\\\\sanus\\\\Scripts\\\\git\\\\tissue-specific-FBA-models\\\\Models\\\\maize_root_model.xml\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "io.write_sbml_model(stem_model,\"C:\\\\Users\\\\sanus\\\\Scripts\\\\git\\\\tissue-specific-FBA-models\\\\Models\\\\maize_stem_model.xml\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "io.write_sbml_model(leaf_model_source,\"C:\\\\Users\\\\sanus\\\\Scripts\\\\git\\\\tissue-specific-FBA-models\\\\Models\\\\maize_sourceleaf_model.xml\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "io.write_sbml_model(growing_leaf_model,\"C:\\\\Users\\\\sanus\\\\Scripts\\\\git\\\\tissue-specific-FBA-models\\\\Models\\\\maize_sinkleaf_model.xml\")"
   ]
  },
  {
   "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
}