4. Building a Model

4.1. Building with Metabolite and Reaction

Building a model starts with creating metabolites and reactions. Let’s consider the reaction OAA(abcd) + AcCoA(ef) \(\rightarrow\) Cit(dcbfea) in the toy model as an example. To create the reactants, we can use the following code:

from freeflux import Metabolite, Reaction, Model

oaa = Metabolite('OAA', atoms = ['a', 'b', 'c', 'd'])
accoa = Metabolite('AcCoA', atoms = ['d', 'f'])
cit = Metabolite('Cit', atoms = list('dcbfea'))
print(oaa, accoa, cit)

Metabolite OAA(a,b,c,d) Metabolite AcCoA(d,f) Metabolite Cit(d,c,b,f,e,a)

Note: Rotationally symmetric metabolites usually have equivalents in biochemical reactions, e.g., succinate and fumarate in the TCA cycle. In this case, the “atoms” argument in Metabolite constructor needs to be assigned as “abcd,dcba”.

Then, we can build a reaction that consumes and produces these metabolites:

v1 = Reaction('v1', reversible = False)
v1.add_substrates([oaa, accoa], stoichiometry = [1, 1])
v1.add_products(cit, stoichiometry = 1)
print(v1)

Reaction v1: AcCoA+OAA->Cit

Finally, we add the reaction to the model:

demo = Model('demo')
demo.add_reactions([v1])
print(demo)

Model demo (3 metabolites, 1)
Reaction v1: AcCoA+OAA->Cit

The process can be repeated until all reactions are involved.

To modify the mode, one can use remove_substrates and remove_products to delete reactants from a reaction, and remove_reactions to delate reactions from a model.

4.2. Reading from File

To input a large metabolic network, it is convenient to load from a file. The following formats are supported: tab-separated values (.tsv) or Excel spreadsheet (.xlsx). The file should have the following format:

#reaction_ID	reactant_IDs(atom)	product_IDs(atom)	reversibility
v1	OAA(abcd)+AcCoA(ef)	Cit(dcbfea)	0
v2	Cit(abcdef)	AKG(abcde)+CO2(f)	0
v3	AKG(abcde)	Glu(abcde)	0
v4	AKG(abcde)	Suc(bcde)+CO2(a)	0
v5	Suc(abcd,dcba)	Fum(abcd,dcba)	0
v6	Fum(abcd,dcba)	OAA(abcd)	1
v7	Asp(abcd)	OAA(abcd)	0

Note: 1. “#” is required in the header; 2. Metabolite name could include but must not start with digits; 3. Reactions with end metabolite, i.e., initial substrates and final products of the network (for example, substrate glucose and biomass) should be irreversible.

To load the model from the file, use the read_from_file method:

MODEL_FILE = 'path/to/reactions.tsv'

demo.read_from_file(MODEL_FILE)
print(demo)

Model demo (9 metabolites, 7)
Reaction v1: AcCoA+OAA->Cit
Reaction v2: Cit->AKG+CO2
Reaction v3: AKG->Glu
Reaction v4: AKG->CO2+Suc
Reaction v5: Suc->Fum
Reaction v6: Fum<->OAA
Reaction v7: Asp->OAA

The model file can be found here.

The metabolites_info and reactions_info attributes can be used to access information about the metabolites and reactions in the model:

demo.metabolites_info

{'OAA': [Metabolite OAA(abcd)],
 'AcCoA': [Metabolite AcCoA(ef)],
 'Cit': [Metabolite Cit(dcbfea), Metabolite Cit(abcdef)],
 'AKG': [Metabolite AKG(abcde)],
 'CO2': [Metabolite CO2(f), Metabolite CO2(a)],
 'Glu': [Metabolite Glu(abcde)],
 'Suc': [Metabolite Suc(bcde), Metabolite Suc(abcd,dcba)],
 'Fum': [Metabolite Fum(abcd,dcba)],
 'Asp': [Metabolite Asp(abcd)]}

demo.reactions_info

OrderedDict([('v1', Reaction v1: AcCoA+OAA->Cit),
             ('v2', Reaction v2: Cit->AKG+CO2),
             ('v3', Reaction v3: AKG->Glu),
             ('v4', Reaction v4: AKG->CO2+Suc),
             ('v5', Reaction v5: Suc->Fum),
             ('v6', Reaction v6: Fum<->OAA),
             ('v7', Reaction v7: Asp->OAA)])

To obtain the stoichiometric matrix for the net reactions, use the get_net_stoichiometric_matrix method:

demo.get_net_stoichiometric_matrix()

	v1	v2	v3	v4	v5	v6	v7
AKG	0.0	1.0	-1.0	-1.0	0.0	0.0	0.0
Cit	1.0	-1.0	0.0	0.0	0.0	0.0	0.0
Fum	0.0	0.0	0.0	0.0	1.0	-1.0	0.0
OAA	-1.0	0.0	0.0	0.0	0.0	1.0	1.0
Suc	0.0	0.0	0.0	1.0	-1.0	0.0	0.0

Similarly, the stoichiometric matrix for the total reactions can be obtained using the get_total_stoichiometric_matrix method:

demo.get_total_stoichiometric_matrix()

	v1	v2	v3	v4	v5	v6_f	v6_b	v7
AKG	0.0	1.0	-1.0	-1.0	0.0	0.0	-0.0	0.0
Cit	1.0	-1.0	0.0	0.0	0.0	0.0	-0.0	0.0
Fum	0.0	0.0	0.0	0.0	1.0	-1.0	1.0	0.0
OAA	-1.0	0.0	0.0	0.0	0.0	1.0	-1.0	1.0
Suc	0.0	0.0	0.0	1.0	-1.0	0.0	-0.0	0.0

4.3. Metabolic Network Decomposition

In order to establish the mapping between fluxes and labeling patterns of metabolites, or more specifically, the elementary metabolite units (EMU), a metabolic network carrying atom transfer information needs to be decomposed. This can be achieved using the adjacency matrix enabled method that has been previously proposed and is implemented by the decompose_network function. To decompose the network, one can use the following code:

emu_mats = demo.decompose_network(
    ini_emus = {'Glu': ['12345']}   # define initial EMU Glu_12345 for decomposition
)

The function returns a dictionary of the decomposed EMU networks with the size as the key.

Note: The decompose_network can run with parallel jobs for multiple target EMUs by specifing the “n_jobs” argument.

For example, to obtain the EMU adjacency matrix for EMUs of size 1, the following code can be used:

emu_mats[1]

	EMU Fum_2	EMU OAA_2	EMU OAA_3
EMU Fum_2	0	1.0*v6_f	1.0*v6_f
EMU OAA_2	0.5*v5 + 0.5*v6_b	0	0
EMU OAA_3	0.5*v6_b	0	0
EMU AcCoA_2	0.5*v5	0	0
EMU Asp_2	0	1.0*v7	0
EMU Asp_3	0	0	1.0*v7

This EMU adjacency matrix shows the transformation of size-1 EMUs bridged by reactions. For instance, the EMU OAA_2 comes from precursor EMU Fum_2 through the forward reaction v6_f and the substrate EMU Asp_2 through the reaction v7. Through the decomposition, the labeling pattern of any target EMU(s) can be traced back to the initial labeled or unlabeled substrate(s).