Building model

This is the code that is used to build the 4-phase integrated guard cell and mesophyll cell model from the original core plant model

Set up

import cobra
import pandas as pd
from mmon_gcm.buildingediting import (
    add_linkers,
    add_metabolite,
    add_metabolites_to_reaction_multi,
    add_reaction,
    check_number_of_models,
    set_bounds_multi,
    split_model,
)

model = cobra.io.read_sbml_model("../models/PlantCoreMetabolism_v1_2_2_3fbc.xml")

No objective coefficients in model. Unclear what should be optimized

Add charge to model

charges = pd.read_csv("../inputs/charges.csv")
charges.set_index("Metabolite")
for index, row in charges.iterrows():
    model.metabolites.get_by_id(row[1]).charge = row[0]

Changes to Both Cells

K+

for metabolite in model.metabolites:
    if "KI_" in metabolite.id:
        metabolite.id = metabolite.id[:1] + metabolite.id[-2:]
        metabolite.name = metabolite.id
model.reactions.K_rev_vc.name = "K_cv"
model.reactions.K_rev_vc.id = "K_cv"
model.metabolites.K_e.charge = 1

Cl-

add_metabolite(model, "Cl_e", "e", multi=False)
add_metabolite(model, "Cl_c", "c", multi=False)
add_metabolite(model, "Cl_v", "v", multi=False)
model.metabolites.Cl_e.notes = {"INCHI": ["InChI=1S/ClH/h1H/p-1"], "SMILES": ["Cl-"]}
model.metabolites.Cl_e.charge = -1
model.metabolites.Cl_c.notes = model.metabolites.Cl_e.notes.copy()
model.metabolites.Cl_c.charge = -1
model.metabolites.Cl_v.notes = model.metabolites.Cl_e.notes.copy()
model.metabolites.Cl_v.charge = -1

add_reaction(model, "Cl_tx", multi=False)
model.reactions.Cl_tx.add_metabolites({"Cl_e": 1})
model.reactions.Cl_tx.lower_bound = -1000

CIT

for metabolite in model.metabolites:
    if "CIT_" in metabolite.id and "aCIT" not in metabolite.id:
        metabolite.charge = -3

model.metabolites.aCIT_v.charge = -2

MAL

for metabolite in model.metabolites:
    if "MAL_" in metabolite.id and "aMAL" not in metabolite.id:
        metabolite.charge = -2

model.metabolites.aMAL_v.charge = -1

add_metabolite(model, "MAL_e", "e", multi=False)
add_metabolite(model, "aMAL_e", "e", multi=False)

Metabolite identifier	aMAL_e
Name	aMAL_e
Memory address	0x7e0687d07f40
Formula	None
Compartment	e
In 0 reaction(s)

model.metabolites.MAL_e.notes = model.metabolites.MAL_v.notes.copy()
model.metabolites.aMAL_e.notes = model.metabolites.aMAL_v.notes.copy()
model.metabolites.MAL_e.charge = model.metabolites.MAL_v.charge
model.metabolites.aMAL_e.charge = model.metabolites.aMAL_v.charge

add_reaction(model, "MAL_tx", multi=False)
model.reactions.MAL_tx.add_metabolites(
    {
        "MAL_e": 0.7,
        "aMAL_e": 0.3,
    }
)
model.reactions.MAL_tx.lower_bound = -1000

MAL Anion Export Channel (Cl/Mal)

add_reaction(model, "MAL_ce", multi=False)
model.reactions.MAL_ce.name = "MAL R/S-Type Anion Channel"
model.reactions.MAL_ce.add_metabolites(
    {"MAL_c": -1, "MAL_e": 0.7, "aMAL_e": 0.3, "PROTON_e": -0.3}
)

NADPHox

# change names to make consistent with the rest of the model
for reaction in model.reactions:
    if "NADPHox" in reaction.id:
        reaction.name = reaction.id[:7] + "_" + reaction.id[7:]
        reaction.id = reaction.id[:7] + "_" + reaction.id[7:]

Remove reactions only active in germinating seeds

Arabidopsis Peroxisomal Citrate Synthase Is Required for Fatty Acid Respiration and Seed Germination Itsara Pracharoenwattana, Johanna E. Cornah, Steven M. Smith The Plant Cell Jul 2005, 17 (7) 2037-2048; DOI: 10.1105/tpc.105.031856

Coordinate expression of transcriptionally regulated isocitrate lyase and malate synthase genes in Brassica napus L. L Comai, R A Dietrich, D J Maslyar, C S Baden, J J Harada The Plant Cell Mar 1989, 1 (3) 293-300; DOI: 10.1105/tpc.1.3.293

model.reactions.MALSYN_RXN_x.bounds = 0, 0
model.reactions.CITSYN_RXN_x.bounds = 0, 0

Setting disproportionating enzyme irreversible in the direction of maltose degradation

Critchley, J.H., Zeeman, S.C., Takaha, T., Smith, A.M. and Smith, S.M. (2001), A critical role for disproportionating enzyme in starch breakdown is revealed by a knock-out mutation in Arabidopsis. The Plant Journal, 26: 89-100. https://doi.org/10.1046/j.1365-313x.2001.01012.x

model.reactions.MALTODEG_RXN_c.bounds = 0, cobra.Configuration().upper_bound

Split Model into Mesophyll (me) and Guard Cell (gc)

Usesplit_model to duplicate the model and add the tags “me” and “gc” to each one

split_model(model, ["me", "gc"])

/home/maurice/miniconda3/envs/mmon-gcm/lib/python3.9/site-packages/cobra/core/group.py:147: UserWarning: need to pass in a list
  warn("need to pass in a list")

Adding GC-specific transporters

model.reactions.PROTON_ATPase_c_gc

Reaction identifier	PROTON_ATPase_c_gc
Name	PROTON_ATPase_c_gc
Memory address	0x7e0686ad35b0
Stoichiometry	0.65 ATP_c_gc + 0.45 PROTON_c_gc + WATER_c_gc + 0.35 aATP_c_gc --> 0.5 ADP_c_gc + PROTON_e_gc + 0.7 Pi_c_gc + 0.5 aADP_c_gc + 0.3 aPi_c_gc 0.65 ATP_gc + 0.45 PROTON_gc + WATER_gc + 0.35 aATP[c]_gc --> 0.5 ADP_gc + PROTON_gc + 0.7 Pi[c]_gc + 0.5 aADP[c]_gc + 0.3 aPi[c]_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Inward-rectifying K+ Channel

Remove dependence for protons

model.reactions.K_ec_gc.name = "Apoplastic Inward-Rectifying K+ Channel"
model.reactions.K_ec_gc.add_metabolites({"PROTON_e_gc": 1, "PROTON_c_gc": -1})

model.reactions.K_ec_gc

Reaction identifier	K_ec_gc
Name	Apoplastic Inward-Rectifying K+ Channel
Memory address	0x7e0686b12bb0
Stoichiometry	K_e_gc --> K_c_gc K_e_gc --> K_c_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Outward-rectifying K+ Channel

add_reaction(model, "K_ce_gc", multi=False)
model.reactions.K_ce_gc.name = "Apoplastic Outward-Rectifying K+ Channel"
model.reactions.K_ce_gc.add_metabolites({"K_c_gc": -1, "K_e_gc": 1})

H+-coupled K+ Symport

add_reaction(model, "K_PROTON_ec_gc", multi=False)
model.reactions.K_PROTON_ec_gc.name = "H+-Coupled K+ Symport"
model.reactions.K_PROTON_ec_gc.add_metabolites(
    {"K_e_gc": -1, "K_c_gc": 1, "PROTON_e_gc": -1, "PROTON_c_gc": 1}
)
model.reactions.K_PROTON_ec_gc

Reaction identifier	K_PROTON_ec_gc
Name	H+-Coupled K+ Symport
Memory address	0x7e0687d43760
Stoichiometry	K_e_gc + PROTON_e_gc --> K_c_gc + PROTON_c_gc K_e_gc + PROTON_gc --> K_c_gc + PROTON_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Cl- Import Channel

add_reaction(model, "Cl_ec_gc", multi=False)
model.reactions.Cl_ec_gc.name = "Cl Apoplastic Import Channel"
model.reactions.Cl_ec_gc.add_metabolites({"Cl_e_gc": -1, "Cl_c_gc": 1})

H+-Coupled Cl- Symport

add_reaction(model, "Cl_PROTON_ec_gc", multi=False)
model.reactions.Cl_PROTON_ec_gc.name = "H+-Coupled Cl- Symport"
model.reactions.Cl_PROTON_ec_gc.add_metabolites(
    {"Cl_e_gc": -2, "Cl_c_gc": 2, "PROTON_e_gc": -1, "PROTON_c_gc": 1}
)

Cl- Anion Export Channel (Cl/Mal)

add_reaction(model, "Cl_ce_gc", multi=False)
model.reactions.Cl_ce_gc.name = "Cl- R/S-Type Anion Channel"
model.reactions.Cl_ce_gc.add_metabolites({"Cl_c_gc": -1, "Cl_e_gc": 1})

ATPase Malate Importer

add_reaction(model, "MAL_ATPASE_ec_gc", multi=False)
model.reactions.MAL_ATPASE_ec_gc.add_metabolites(
    {
        "MAL_e_gc": -0.7,
        "aMAL_e_gc": -0.3,
        "MAL_c_gc": 1,
        "ATP_c_gc": -0.65,
        "aATP_c_gc": -0.35,
        "WATER_c_gc": -1,
        "ADP_c_gc": 1,
        "Pi_c_gc": 0.7,
        "aPi_c_gc": 0.3,
        "PROTON_c_gc": +0.85,
    }
)

Glucose Apoplastic Symport Channel

model.reactions.GLC_ec_gc.id = "GLC_PROTON_ec_gc"
model.reactions.GLC_PROTON_ec_gc.name = "H+-Coupled Glucose Symport"
model.reactions.GLC_PROTON_ec_gc

Reaction identifier	GLC_PROTON_ec_gc
Name	H+-Coupled Glucose Symport
Memory address	0x7e0686c9c8e0
Stoichiometry	GLC_e_gc + PROTON_e_gc --> GLC_c_gc + PROTON_c_gc GLC_gc + PROTON_gc --> GLC_gc + PROTON_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Glucose Exporter

add_reaction(model, "GLC_ce_gc", multi=False)
model.reactions.GLC_ce_gc.name = "Glucose Apoplastic Exporter"
model.reactions.GLC_ce_gc.add_metabolites({"GLC_c_gc": -1, "GLC_e_gc": 1})

Sucrose Apoplastic Symport Channel

# Renaming Sucrose_tx
model.reactions.Sucrose_tx_gc.id = "SUCROSE_tx_gc"

model.reactions.Sucrose_ec_gc.id = "SUCROSE_PROTON_ec_gc"
model.reactions.SUCROSE_PROTON_ec_gc.name = "H+-Coupled Sucrose Symport"
model.reactions.SUCROSE_PROTON_ec_gc.upper_bound = 1000
model.reactions.SUCROSE_PROTON_ec_gc

Reaction identifier	SUCROSE_PROTON_ec_gc
Name	H+-Coupled Sucrose Symport
Memory address	0x7e06869a2cd0
Stoichiometry	PROTON_e_gc + SUCROSE_e_gc --> PROTON_c_gc + SUCROSE_c_gc PROTON_gc + SUCROSE_gc --> PROTON_gc + SUCROSE_gc
GPR
Lower bound	0.0
Upper bound	1000

Sucrose Exporter

add_reaction(model, "SUCROSE_ce_gc", multi=False)
model.reactions.SUCROSE_ce_gc.name = "Sucrose Apoplastic Exporter"
model.reactions.SUCROSE_ce_gc.add_metabolites({"SUCROSE_c_gc": -1, "SUCROSE_e_gc": 1})

Fructose Apoplastic Symport Channel

add_metabolite(model, "FRU_e_gc", "e", multi=False)
model.metabolites.FRU_e_gc.notes = model.metabolites.FRU_c_gc.notes.copy()
model.metabolites.FRU_e_gc.charge = 0

add_reaction(model, "FRU_PROTON_ec_gc", multi=False)
model.reactions.FRU_PROTON_ec_gc.name = "H+-Coupled Fructose Symport"
model.reactions.FRU_PROTON_ec_gc.add_metabolites(
    {"FRU_e_gc": -1, "FRU_c_gc": 1, "PROTON_e_gc": -1, "PROTON_c_gc": 1}
)

Fructose Exporter

add_reaction(model, "FRU_ce_gc", multi=False)
model.reactions.FRU_ce_gc.name = "Fructose Apoplastic Exporter"
model.reactions.FRU_ce_gc.add_metabolites({"FRU_c_gc": -1, "FRU_e_gc": 1})

Cell-Wall Invertase

add_reaction(model, "cwINV_gc", multi=False)

Reaction identifier	cwINV_gc
Name	cwINV_gc
Memory address	0x7e0687d43a60
Stoichiometry	--> -->
GPR
Lower bound	0.0
Upper bound	1000.0

model.reactions.cwINV_gc.add_metabolites(
    {"SUCROSE_e_gc": -1, "WATER_e_gc": -1, "FRU_e_gc": 1, "GLC_e_gc": 1}
)

K+ Import Channel

(TPK-Type and Fast-vacuolar)

model.reactions.K_cv_gc.name = "K+ Tonoplastic Import Channel"
model.reactions.K_cv_gc

Reaction identifier	K_cv_gc
Name	K+ Tonoplastic Import Channel
Memory address	0x7e0686aa7dc0
Stoichiometry	K_c_gc + PROTON_v_gc --> K_v_gc + PROTON_c_gc K_c_gc + PROTON_gc --> K_v_gc + PROTON_gc
GPR
Lower bound	0.0
Upper bound	1000.0

NHX-Type K+/H+ Antiport

add_reaction(model, "K_PROTON_cv_gc", multi=False)
model.reactions.K_PROTON_cv_gc.name = "NHX-Type K+/H+ Antiport"
model.reactions.K_PROTON_cv_gc.add_metabolites(
    {"K_c_gc": -1, "K_v_gc": 1, "PROTON_v_gc": -1, "PROTON_c_gc": 1}
)

K+ Tonoplastic Export Channel

Slow Vacuolar Channel

add_reaction(model, "K_vc_gc", multi=False)
model.reactions.K_vc_gc.name = "K+ Tonoplastic Export Slow Vacuolar Channel"
model.reactions.K_vc_gc.add_metabolites({"K_v_gc": -1, "K_c_gc": 1})

CLC-Type Cl-/H+ Antiport Vacoular Import Channel

add_reaction(model, "Cl_PROTON_cv_gc", multi=False)
model.reactions.Cl_PROTON_cv_gc.name = (
    "CLC-Type Cl-/H+ Antiport Vacoular Import Channel"
)
model.reactions.Cl_PROTON_cv_gc.add_metabolites(
    {
        "Cl_v_gc": 1,
        "Cl_c_gc": -1,
        "PROTON_c_gc": 2,
        "PROTON_v_gc": -2,
    }
)

VCL Cl- Vacuolar Import Channel

add_reaction(model, "Cl_cv_gc", multi=False)
model.reactions.Cl_cv_gc.name = "VCL Cl- Vacuolar Import Channel"
model.reactions.Cl_cv_gc.add_metabolites({"Cl_c_gc": -1, "Cl_v_gc": 1})

Cl- Vacuolar Export Channel

model.reactions.Cl_cv_gc.lower_bound = -1000

MAL Import Channel (AMLT)

model.reactions.MAL_PROTON_vc_gc.id = "MAL_cv_gc"
model.reactions.MAL_cv_gc.name = "VMAL-type MAL Channel (Import)"
model.reactions.MAL_cv_gc

Reaction identifier	MAL_cv_gc
Name	VMAL-type MAL Channel (Import)
Memory address	0x7e0686ab9190
Stoichiometry	MAL_c_gc + 0.3 PROTON_v_gc --> 0.7 MAL_v_gc + 0.3 aMAL_v_gc MAL_gc + 0.3 PROTON_gc --> 0.7 MAL_gc + 0.3 aMAL[v]_gc
GPR
Lower bound	0.0
Upper bound	1000.0

MAL Export Channel (AMLT)

add_reaction(model, "MAL_vc_gc", multi=False)
model.reactions.MAL_vc_gc.name = "VMAL-type MAL Channel (Export)"
model.reactions.MAL_vc_gc.add_metabolites(
    {
        "MAL_v_gc": -0.7,
        "aMAL_v_gc": -0.3,
        "MAL_c_gc": 1,
        "PROTON_c_gc": 0.3,
    }
)

MAL Export Symporter

Already in Model

model.reactions.MAL_PROTON_rev_vc_gc.name = "MAL Export Symporter"
model.reactions.MAL_PROTON_rev_vc_gc

Reaction identifier	MAL_PROTON_rev_vc_gc
Name	MAL Export Symporter
Memory address	0x7e068696ca90
Stoichiometry	0.7 MAL_v_gc + 1.7 PROTON_v_gc + 0.3 aMAL_v_gc --> MAL_c_gc + 2.0 PROTON_c_gc 0.7 MAL_gc + 1.7 PROTON_gc + 0.3 aMAL[v]_gc --> MAL_gc + 2.0 PROTON_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Fructose Import Antiporter

Previously in the model

model.reactions.FRU_PROTON_rev_vc_gc.id = "FRU_PROTON_rev_cv_gc"
model.reactions.FRU_PROTON_rev_cv_gc.name = "Fructose Tonoplastic Import Antiporter"
model.reactions.FRU_PROTON_rev_cv_gc

Reaction identifier	FRU_PROTON_rev_cv_gc
Name	Fructose Tonoplastic Import Antiporter
Memory address	0x7e0686ea8fa0
Stoichiometry	FRU_c_gc + PROTON_v_gc --> FRU_v_gc + PROTON_c_gc FRU[c]_gc + PROTON_gc --> FRU[v]_gc + PROTON_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Fructose Import Symporter

add_reaction(model, "FRU_PROTON_cv_gc", multi=False)
model.reactions.FRU_PROTON_cv_gc.name = "Fructose Tonoplastic Import Symporter"
model.reactions.FRU_PROTON_cv_gc.add_metabolites(
    {"PROTON_c_gc": -1, "PROTON_v_gc": 1, "FRU_c_gc": -1, "FRU_v_gc": 1}
)

Fructose Export

model.reactions.FRU_vc_gc

Reaction identifier	FRU_vc_gc
Name	FRU_vc_gc
Memory address	0x7e0686b245e0
Stoichiometry	FRU_v_gc --> FRU_c_gc FRU[v]_gc --> FRU[c]_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Glucose Tonoplastic Antiporter

model.reactions.GLC_PROTON_rev_vc_gc.id = "GLC_PROTON_rev_cv_gc"
model.reactions.GLC_PROTON_rev_cv_gc.name = "Glucose Tonoplastic Import Antiporter"
model.reactions.GLC_PROTON_rev_cv_gc

Reaction identifier	GLC_PROTON_rev_cv_gc
Name	Glucose Tonoplastic Import Antiporter
Memory address	0x7e0686d1de50
Stoichiometry	GLC_c_gc + PROTON_v_gc --> GLC_v_gc + PROTON_c_gc GLC_gc + PROTON_gc --> GLC_gc + PROTON_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Glucose Tonoplastic Symporter

add_reaction(model, "GLC_PROTON_cv_gc", multi=False)
model.reactions.GLC_PROTON_cv_gc.name = "Glucose Tonoplastic Import Symporter"
model.reactions.GLC_PROTON_cv_gc.add_metabolites(
    {"PROTON_c_gc": -1, "PROTON_v_gc": 1, "GLC_c_gc": -1, "GLC_v_gc": 1}
)

Glucose Export Channel

model.reactions.GLC_vc_gc

Reaction identifier	GLC_vc_gc
Name	GLC_vc_gc
Memory address	0x7e0686a9df70
Stoichiometry	GLC_v_gc --> GLC_c_gc GLC_gc --> GLC_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Sucrose Tonoplastic Importer

add_reaction(model, "SUCROSE_cv_gc", multi=False)
model.reactions.SUCROSE_cv_gc.name = "Sucrose Free Tonoplastic Import"
model.reactions.SUCROSE_cv_gc.add_metabolites(
    {
        "SUCROSE_c_gc": -1,
        "SUCROSE_v_gc": 1,
    }
)

Sucrose Tonoplastic Import Antiporter

Already in the model

model.reactions.SUCROSE_PROTON_vc_gc.id = "SUCROSE_PROTON_cv_gc"
model.reactions.SUCROSE_PROTON_cv_gc.name = "Sucrose Tonoplastic Import Antiporter"
model.reactions.SUCROSE_PROTON_cv_gc

Reaction identifier	SUCROSE_PROTON_cv_gc
Name	Sucrose Tonoplastic Import Antiporter
Memory address	0x7e0686f28d30
Stoichiometry	PROTON_v_gc + SUCROSE_c_gc --> PROTON_c_gc + SUCROSE_v_gc PROTON_gc + SUCROSE_gc --> PROTON_gc + SUCROSE_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Sucrose Tonoplastic Export Symporter

model.reactions.SUCROSE_PROTON_rev_vc_gc.id = "SUCROSE_PROTON_vc_gc"
model.reactions.SUCROSE_PROTON_vc_gc.name = "Sucrose Tonoplastic Export Symporter"

Tonoplastic PPase

model.reactions.PROTON_PPi_rev_vc_gc

Reaction identifier	PROTON_PPi_rev_vc_gc
Name	PROTON_PPi_rev_vc_gc
Memory address	0x7e06869e5670
Stoichiometry	0.65 PPI_c_gc + 0.25 PROTON_c_gc + WATER_c_gc + 0.35 aPPI_c_gc --> PROTON_v_gc + 1.4 Pi_c_gc + 0.6 aPi_c_gc 0.65 PPI_gc + 0.25 PROTON_gc + WATER_gc + 0.35 aPPI[c]_gc --> PROTON_gc + 1.4 Pi[c]_gc + 0.6 aPi[c]_gc
GPR
Lower bound	0.0
Upper bound	1000.0

Tonoplastic H+ ATPase

model.reactions.PROTONATP_rev_vc_gc

Reaction identifier	PROTONATP_rev_vc_gc
Name	PROTONATP_rev_vc_gc
Memory address	0x7e0686febb20
Stoichiometry	0.65 ATP_c_gc + 1.45 PROTON_c_gc + WATER_c_gc + 0.35 aATP_c_gc --> 0.5 ADP_c_gc + 2.0 PROTON_v_gc + 0.7 Pi_c_gc + 0.5 aADP_c_gc + 0.3 aPi_c_gc 0.65 ATP_gc + 1.45 PROTON_gc + WATER_gc + 0.35 aATP[c]_gc --> 0.5 ADP_gc + 2.0 PROTON_gc + 0.7 Pi[c]_gc + 0.5 aADP[c]_gc + 0.3 aPi[c]_gc
GPR
Lower bound	0.0
Upper bound	1000.0

add_reaction(model, "PROTON_ec_gc", multi=False)
model.reactions.PROTON_ec_gc.name = "Plasma membrane proton leakage"
model.reactions.PROTON_ec_gc.add_metabolites(
    {
        "PROTON_e_gc": -1,
        "PROTON_c_gc": 1,
    }
)

add_reaction(model, "PROTON_vc_gc", multi=False)
model.reactions.PROTON_vc_gc.name = "Tonoplast membrane proton leakage"
model.reactions.PROTON_vc_gc.add_metabolites(
    {
        "PROTON_v_gc": -1,
        "PROTON_c_gc": 1,
    }
)

Constrain Channels in Guard Cell

for channel in ["Cl_ec", "K_cv", "FRU_PROTON_cv", "GLC_PROTON_cv", "SUCROSE_cv"]:
    model.reactions.get_by_id(channel + "_gc").bounds = 0, 0

Adding apoplast and exchanges

In this section we add an apoplastic compartment and set the influx of osmolites to be from this compartment rather than as boundary reactions.

Removing tx_me

Here we remove all the boundary reactions from the mesophyll cell except for: - Maintenance reactions - Photon influx - Phloem output - Gas transfer (O2 and CO2)

tx_me = []
for reaction in model.reactions:
    if "_tx_me" in reaction.id:
        tx_me.append(reaction)
keep = [
    "Photon_tx_me",
    "Phloem_output_tx_me",
    "NADPHox_c_tx_me",
    "NADPHox_m_tx_me",
    "NADPHox_p_tx_me",
    "O2_tx_me",
    "CO2_tx_me",
    "ATPase_tx_me",
]
for reaction in keep:
    tx_me.remove(model.reactions.get_by_id(reaction))

model.remove_reactions(tx_me)

Removing tx_gc

See above for removing tx_me reactions, only difference is in this case we don’t have a phloem output reaction

tx_gc = []
for reaction in model.reactions:
    if "_tx_gc" in reaction.id:
        tx_gc.append(reaction)
keep = [
    "Photon_tx_gc",
    "NADPHox_c_tx_gc",
    "NADPHox_m_tx_gc",
    "NADPHox_p_tx_gc",
    "O2_tx_gc",
    "CO2_tx_gc",
    "ATPase_tx_gc",
]
for reaction in keep:
    tx_gc.remove(model.reactions.get_by_id(reaction))

model.remove_reactions(tx_gc)

Adding tx to apo and adding exchanges

Now for the species that are present in the extracellular compartment of the guard cell we: - Add the species to the apoplastic compartment with the tag “a” - Add a boundary reaction into this compartment for the species - Add a transfer reaction to allow them to be transferred between the guard cell extracellular compartment and the apoplast - Add a transfer reaction to allow them to be transferred between the mesophyll cell extracellular compartment and the apoplast

All these actions are free and reversible, as we are assuming that really the extracellular space and the apoplast are all the same. The apoplastic compartment just gives us an easy way to track influx/efflux of metabolites as well as osmolite levels in that compartment.

# fructose wasn't in the extracellular compartment of mesophyll before
add_metabolite(model, "FRU_e_me", "e", multi=False)

Metabolite identifier	FRU_e_me
Name	FRU_e_me
Memory address	0x7e0687e36250
Formula	None
Compartment	e
In 0 reaction(s)

add_reaction(model, "Sucrose_ce_me", multi=False)  # sucrose transport out of the mesophyll cell wasn't allowed before
model.reactions.Sucrose_ce_me.add_metabolites({"SUCROSE_c_me": -1, "SUCROSE_e_me": 1})

ea_reactions = []
for metabolite in model.metabolites:
    if "_e_gc" in metabolite.id:
        ea_reactions.append(metabolite)

print(ea_reactions)

# don't want these to be transferred in the apoplast, they already have individual boundary reactions in each cell type
remove = ["Photon", "OXYGEN_MOLECULE", "CARBON_DIOXIDE", "PROTON"]
for name in remove:
    ea_reactions.remove(model.metabolites.get_by_id(name + "_e_gc"))

for metabolite in ea_reactions:

    add_metabolite(model, metabolite.id[:-4] + "a", "a", multi=False)
    model.metabolites.get_by_id(metabolite.id[:-4] + "a").charge = metabolite.charge
    model.metabolites.get_by_id(metabolite.id[:-4] + "a").notes = metabolite.notes

    add_reaction(model, metabolite.id[:-4] + "a_tx", multi=False)
    model.reactions.get_by_id(metabolite.id[:-4] + "a_tx").add_metabolites(
        {metabolite.id[:-4] + "a": 1}
    )
    model.reactions.get_by_id(metabolite.id[:-4] + "a_tx").lower_bound = -1000

    add_reaction(model, metabolite.id[:-4] + "ae_gc", multi=False)
    model.reactions.get_by_id(metabolite.id[:-4] + "ae_gc").add_metabolites(
        {metabolite: 1, metabolite.id[:-4] + "a": -1}
    )
    model.reactions.get_by_id(metabolite.id[:-4] + "ae_gc").lower_bound = -1000

    add_reaction(model, metabolite.id[:-4] + "ae_me", multi=False)
    model.reactions.get_by_id(metabolite.id[:-4] + "ae_me").add_metabolites(
        {metabolite.id[:-4] + "e_me": 1, metabolite.id[:-4] + "a": -1}
    )
    model.reactions.get_by_id(metabolite.id[:-4] + "ae_me").lower_bound = -1000

[<Metabolite NITRATE_e_gc at 0x7e06872203a0>, <Metabolite SUCROSE_e_gc at 0x7e0687220b20>, <Metabolite WATER_e_gc at 0x7e0687236be0>, <Metabolite MGII_e_gc at 0x7e06871c25e0>, <Metabolite GLC_e_gc at 0x7e06871d5b20>, <Metabolite PROTON_e_gc at 0x7e06871dc8e0>, <Metabolite CARBON_DIOXIDE_e_gc at 0x7e068719acd0>, <Metabolite CAII_e_gc at 0x7e06871bd910>, <Metabolite Pi_e_gc at 0x7e06871bdf10>, <Metabolite SULFATE_e_gc at 0x7e0687154250>, <Metabolite OXYGEN_MOLECULE_e_gc at 0x7e0687154fd0>, <Metabolite Photon_e_gc at 0x7e0687119280>, <Metabolite K_e_gc at 0x7e0687125a00>, <Metabolite AMMONIUM_e_gc at 0x7e06870e0580>, <Metabolite Cl_e_gc at 0x7e06870f8b80>, <Metabolite MAL_e_gc at 0x7e06870f8dc0>, <Metabolite aMAL_e_gc at 0x7e06870f8e80>, <Metabolite FRU_e_gc at 0x7e0687c417f0>]

List of all boundary reactions:

for reaction in model.reactions:
    if "_tx" in reaction.id:
        print(reaction.id)

Photon_tx_me
Photon_tx_gc
NADPHox_m_tx_me
NADPHox_m_tx_gc
CO2_tx_me
CO2_tx_gc
O2_tx_me
O2_tx_gc
Phloem_output_tx_me
NADPHox_c_tx_me
NADPHox_c_tx_gc
ATPase_tx_me
ATPase_tx_gc
NADPHox_p_tx_me
NADPHox_p_tx_gc
NITRATE_a_tx
SUCROSE_a_tx
WATER_a_tx
MGII_a_tx
GLC_a_tx
CAII_a_tx
Pi_a_tx
SULFATE_a_tx
K_a_tx
AMMONIUM_a_tx
Cl_a_tx
MAL_a_tx
aMAL_a_tx
FRU_a_tx

Constrain free exchange and starch breakdown

Prevent: - Free import into the apoplast of energy-carrying metabolites like Glucose, malate, sucrose, fructose. - Free exchange of protons. - Alpha-glucosidase starch degradation pathway as it isn’t used

for reaction in [
    "GLC_a_tx",
    "MAL_a_tx",
    "SUCROSE_a_tx",
    "FRU_a_tx",
    "RXN_1826_p_me",
    "RXN_1826_p_gc",
    "MALTODEXGLUCOSID_RXN_p_me",
    "MALTODEXGLUCOSID_RXN_p_gc",
]:
    model.reactions.get_by_id(reaction).bounds = (0, 0)

Quadruple model into 4 phases

split_model(model, range(1, 5))

Add Linker Reactions With Pseudometabolites (Osmolarity and Charge)

Set phase times here

compartments = ["c", "v", "p", "a"]

cells = [
    "gc",
    "me",
]
phase_times = [6.0, 0.5, 11.5, 6.0]

add_linkers(model, "../inputs/osmolytes.csv", compartments, cells, phase_times)

prefix_metabolites = {
    "a": {
        "v_gc": {
            "MAL": 0.7,
            "CIT": 0.5,
        },
        "a": {
            "MAL": 0.7,
        },
        "v_me": {
            "MAL": 0.7,
            "CIT": 0.5,
        },
    },
    "b": {"v_me": {"HIS": 0}, "v_gc": {"HIS": 0}},
}

for prefix, compartments in prefix_metabolites.items():
    for compartment, metabolites in compartments.items():
        for metabolite, ratio in metabolites.items():
            if ratio == 0:
                set_bounds_multi(
                    model, metabolite + "_" + compartment + "_Linker", 0, 0
                )
            elif ratio == 1:
                set_bounds_multi(
                    prefix + model, metabolite + "_" + compartment + "_Linker", 0, 0
                )
            else:
                add_metabolite(
                    model, metabolite + "_" + compartment + "_prefixpseudometabolite"
                )
                for phase in range(check_number_of_models(model)):
                    model.reactions.get_by_id(
                        f"{metabolite}_{compartment}_Linker_" + str(phase + 1)
                    ).add_metabolites(
                        {
                            model.metabolites.get_by_id(
                                f"{metabolite}_{compartment}_prefixpseudometabolite_" +
                                str(phase + 1)
                            ): 1 -
                            ratio
                        }
                    )
                    model.reactions.get_by_id(
                        f"{prefix}{metabolite}_{compartment}_Linker_" + str(phase + 1)
                    ).add_metabolites(
                        {
                            model.metabolites.get_by_id(
                                f"{metabolite}_{compartment}_prefixpseudometabolite_"
                                + str(phase + 1)
                            ): -ratio
                        }
                    )

Add Maintenance Reactions

Remove maintenance phase constraint (no forcing same maintenance for each phase)

number_of_models = check_number_of_models(model)
for cell in ["gc", "me"]:
    add_metabolite(model, "maintenance_ratio_constraint_" + cell, "pseudo")
    #add_metabolite(model, "maintenance_phase_constraint_" + cell, "pseudo")
    add_reaction(model, "maintenance_phase_overall_" + cell, multi="")
    for i in range(1, number_of_models + 1):
        for x in ["c", "m", "p"]:
            reaction = model.reactions.get_by_id(f"NADPHox_{x}_tx_{cell}_" + str(i))
            reaction.add_metabolites(
                {f"maintenance_ratio_constraint_{cell}_" + str(i): -3}
            )
        reaction = model.reactions.get_by_id(f"ATPase_tx_{cell}_" + str(i))
        reaction.add_metabolites(
            {
                f"maintenance_ratio_constraint_{cell}_" + str(i): 1,
                #f"maintenance_phase_constraint_{cell}_" + str(i): 1,
            }
        )
        #model.reactions.get_by_id(f"maintenance_phase_overall_{cell}").add_metabolites(
        #    {f"maintenance_phase_constraint_{cell}_" + str(i): -1}
        #)

Add Phloem_tx Overall

number_of_models = check_number_of_models(model)
day = [2, 3]
night = [1, 4]
add_metabolite(model, "pseudoPhloem_me", "pseudo")
add_metabolite(model, "pseudoPhloem_day_me", "pseudo", multi="")
add_metabolite(model, "pseudoPhloem_night_me", "pseudo", multi="")
add_reaction(model, "Phloem_constraint_day", multi="")
add_reaction(model, "Phloem_constraint_night", multi="")
add_reaction(model, "Phloem_tx_overall", multi="")
for i in range(1, number_of_models + 1):
    length_of_phase = 1 / (
        -model.reactions.get_by_id(f"SUCROSE_v_gc_Linker_{i}").get_coefficient(
            f"SUCROSE_v_gc_{i}"
        )
    )
    model.reactions.get_by_id(f"Phloem_output_tx_me_{i}").add_metabolites(
        {f"pseudoPhloem_me_{i}": 1}
    )
    if i in day:
        model.reactions.get_by_id(f"Phloem_output_tx_me_{i}").add_metabolites(
            {"pseudoPhloem_day_me": 1 * length_of_phase}
        )
        model.reactions.Phloem_constraint_day.add_metabolites(
            {f"pseudoPhloem_me_{i}": -1}
        )
    elif i in night:
        model.reactions.get_by_id(f"Phloem_output_tx_me_{i}").add_metabolites(
            {"pseudoPhloem_night_me": 1 * length_of_phase}
        )
        model.reactions.Phloem_constraint_night.add_metabolites(
            {f"pseudoPhloem_me_{i}": -1}
        )
    else:
        raise ValueError("Make sure all phases are either assigned to day or night")
model.reactions.Phloem_constraint_day.add_metabolites({})
model.reactions.Phloem_tx_overall.add_metabolites(
    {
        "pseudoPhloem_day_me": -3,
        "pseudoPhloem_night_me": -1,
    }
)

Add pseudoreactions for total metabolite counting and FVA

Essentially, here we add a ‘total’ pseudometabolite in the guard cel for each of a list of metabolites. We can use this metabolite to see the sum of the linker reactions for the vacuole and cytoplasm, but most importantly we can perform FVA on this reaction to see if the total transfer of a metabolite matters. Otherwise, we end up in a situation where we can do FVA on the cytoplasmic and vacuolar reactions separately but it could be the case that when cytoplasm goes down the vacuole compensates etc.

"total" added as keyword to getweightings in order to prevent these reactions from contributing to sum of fluxes constraint

metabolites = ["SUCROSE", "GLC", "MAL", "FRU", "K", "Cl", "CIT"]  # "NITRATE", "CIT"]


for metabolite in metabolites:
    if metabolite == "MAL" or metabolite == "CIT":
        add_reaction(model, f"{metabolite}_total_pseudolinker")
        add_metabolite(model, f"{metabolite}_total_pseudometabolite")
        add_metabolites_to_reaction_multi(
            model,
            f"{metabolite}_total_pseudolinker",
            {f"{metabolite}_total_pseudometabolite": -1},
        )
        add_metabolites_to_reaction_multi(
            model,
            f"{metabolite}_c_gc_Linker",
            {f"{metabolite}_total_pseudometabolite": 1},
        )
        add_metabolites_to_reaction_multi(
            model,
            f"{metabolite}_v_gc_Linker",
            {f"{metabolite}_total_pseudometabolite": 1},
        )
        add_metabolites_to_reaction_multi(
            model,
            f"a{metabolite}_v_gc_Linker",
            {f"{metabolite}_total_pseudometabolite": 1},
        )
    else:
        add_reaction(model, f"{metabolite}_total_pseudolinker")
        add_metabolite(model, f"{metabolite}_total_pseudometabolite")
        add_metabolites_to_reaction_multi(
            model,
            f"{metabolite}_total_pseudolinker",
            {f"{metabolite}_total_pseudometabolite": -1},
        )
        add_metabolites_to_reaction_multi(
            model,
            f"{metabolite}_c_gc_Linker",
            {f"{metabolite}_total_pseudometabolite": 1},
        )
        add_metabolites_to_reaction_multi(
            model,
            f"{metabolite}_v_gc_Linker",
            {f"{metabolite}_total_pseudometabolite": 1},
        )

Export Model

cobra.io.write_sbml_model(model, "../models/4_stage_GC.xml")
cobra.io.save_json_model(model, "../models/4_stage_GC.json")