6. Enzyme Modules
An “enzyme module” is defined as a mechanistic description of a reaction consisting of mass action rate laws for all known reaction steps [DZK+16]. In MASSpy, enzyme modules are represented by the EnzymeModule
object.
To demonstrate the utility of an EnzymeModule
object and how it aids in constructing mechanistic models of enzyme behavior, an EnzymeModule
of hexokinase\(^{1, 2}\) is constructed and then merged with a model of glycolysis\(^{3}\) for verification.
6.1. Constructing Enzyme Modules
In order to construct the EnzymeModule
of hexokinase, the following information is necessary:
The enzyme is a monomer.
The enzyme binding of substrates follows a random sequential mechanism.
The enzyme experiences product inhibtion and is competitively inhibited by 23DPG when complexed with D-glucose.
Total HEX1 Concentration\(^2\): \(\text{[HEX1]}_{total} = 24 nM = 0.000024 mM\).
[1]:
from operator import attrgetter
from mass import MassMetabolite
from mass.enzyme_modules import EnzymeModule
from mass.example_data import create_example_model
# Load the glycolysis and hemoglobin models, then merge them
glycolysis = create_example_model("Glycolysis")
hemoglobin = create_example_model("Hemoglobin")
glyc_hb = glycolysis.merge(hemoglobin, inplace=False)
Set parameter Username
The EnzymeModule
is a subclass of the MassModel
, meaning that it inherits the methods and behaviors of the MassModel
object. Like a MassModel
, an EnzymeModule
object requires a unique identifier in order to be created. Optionally, the name
and subsystem
attributes are set during initialization.
[2]:
HEX1 = EnzymeModule("HEX1", name="Hexokinase (D-glucose:ATP)",
subsystem="Glycolysis")
6.1.1. Defining the enzyme ligands
The ligands that interact with the enzyme (e.g. as the substrates, activators, and inhibitors) are created as MassMetabolite
objects and added to the model.
[3]:
glc__D_c = MassMetabolite(
"glc__D_c",
name="D-Glucose",
formula="C6H12O6",
charge=0,
compartment="c")
g6p_c = MassMetabolite(
"g6p_c",
name="D-Glucose 6-phosphate",
formula="C6H11O9P",
charge=-2,
compartment="c")
atp_c = MassMetabolite(
"atp_c",
name="ATP",
formula="C10H12N5O13P3",
charge=-4,
compartment="c")
adp_c = MassMetabolite(
"adp_c",
name="ADP",
formula="C10H12N5O10P2",
charge=-3,
compartment="c")
_23dpg_c = MassMetabolite(
"_23dpg_c",
name="2,3-Disphospho-D-glycerate",
formula="C3H3O10P2",
charge=-5,
compartment="c")
h_c = MassMetabolite(
"h_c",
name="H+",
formula="H",
charge=1,
compartment="c")
HEX1.add_metabolites([glc__D_c, g6p_c, atp_c, adp_c, _23dpg_c, h_c])
Once added to the EnzymeModule
, ligands can be accessed using the enzyme_module_ligands
attribute.
[4]:
HEX1.enzyme_module_ligands
[4]:
[<MassMetabolite glc__D_c at 0x7fb8d158c820>,
<MassMetabolite g6p_c at 0x7fb8d158c7c0>,
<MassMetabolite atp_c at 0x7fb8d158c850>,
<MassMetabolite adp_c at 0x7fb8d158c8b0>,
<MassMetabolite _23dpg_c at 0x7fb8d158c910>,
<MassMetabolite h_c at 0x7fb8d158c940>]
To keep track of the roles played by various ligands in the module, the enzyme_module_ligands_categorized
attribute is set. The attribute takes a dict
, with categories as keys and relevant MassMetabolite
objects as values. Note that an object can be a part of multiple categories.
[5]:
HEX1.enzyme_module_ligands_categorized = {
"substrates": glc__D_c,
"cofactors": atp_c,
"inhibitors": _23dpg_c,
"products": [adp_c, g6p_c, h_c]}
HEX1.enzyme_module_ligands_categorized
[5]:
[<Group substrates at 0x7fb8d158c6d0>,
<Group cofactors at 0x7fb8d158c760>,
<Group inhibitors at 0x7fb8d158cbe0>,
<Group products at 0x7fb8d158cb80>]
For each category, a cobra.Group
is created containing the relevant objects. Once set, the attribute returns a cobra.DictList
that contains the categorized groups. The groups and their members are printed as follows:
[6]:
for group in HEX1.enzyme_module_ligands_categorized:
print("{0}: {1}".format(
group.id, str(sorted([m.id for m in group.members]))))
substrates: ['glc__D_c']
cofactors: ['atp_c']
inhibitors: ['_23dpg_c']
products: ['adp_c', 'g6p_c', 'h_c']
6.1.2. Defining the enzyme module forms
After adding MassMetabolite
objects of ligands to the model, the various forms of the enzyme must be defined. These forms are represented by EnzymeModuleForm
objects.
The EnzymeModuleForm
object inherits from the MassMetabolite
and is treated like any other metabolite in the model. However, the EnzymeModuleForm
object contains the additional bound_metabolites
attribute to assist in tracking metabolites bound to the enzyme form.
The EnzymeModule.make_enzyme_module_form()
method allows for the creation of an EnzymeModuleForm
object while assigning categories for the EnzymeModuleForm
in the process. Using make_enzyme_module_form()
also adds the species to the module upon creation, accessible via the EnzymeModule.enzyme_module_forms
attribute.
[7]:
hex1_c = HEX1.make_enzyme_module_form(
"hex1_c",
name="automatic",
categories="Active",
compartment="c")
hex1_A_c = HEX1.make_enzyme_module_form(
"hex1_A_c", # A stands complexted with ATP
name="automatic",
categories="Active",
bound_metabolites={atp_c: 1},
compartment="c")
hex1_G_c = HEX1.make_enzyme_module_form(
"hex1_G_c", # G stands for complexed with Glucose
name="automatic",
categories="Active",
bound_metabolites={glc__D_c: 1},
compartment="c")
hex1_AG_c = HEX1.make_enzyme_module_form(
"hex1_AG_c",
name="automatic",
categories="Active",
bound_metabolites={glc__D_c: 1, atp_c: 1},
compartment="c")
hex1_G_CI_c = HEX1.make_enzyme_module_form(
"hex1_G_CI_c", # CI stands for competitive inhibition
name="automatic",
categories="Inhibited",
bound_metabolites={glc__D_c: 1, _23dpg_c: 1},
compartment="c")
hex1_A_PI_c = HEX1.make_enzyme_module_form(
"hex1_A_PI_c", # PI stands for competitive inhibition
name="automatic",
categories="Inhibited",
bound_metabolites={adp_c: 1},
compartment="c")
hex1_G_PI_c = HEX1.make_enzyme_module_form(
"hex1_G_PI_c", # PI stands for competitive inhibition
name="automatic",
categories="Inhibited",
bound_metabolites={g6p_c: 1},
compartment="c")
HEX1.enzyme_module_forms
[7]:
[<EnzymeModuleForm hex1_c at 0x7fb8d15a83a0>,
<EnzymeModuleForm hex1_A_c at 0x7fb8d15a8400>,
<EnzymeModuleForm hex1_G_c at 0x7fb8d15a88e0>,
<EnzymeModuleForm hex1_AG_c at 0x7fb8d15a8580>,
<EnzymeModuleForm hex1_G_CI_c at 0x7fb8d15a88b0>,
<EnzymeModuleForm hex1_A_PI_c at 0x7fb8d15a8880>,
<EnzymeModuleForm hex1_G_PI_c at 0x7fb8d15a8160>]
The bound_metabolites
attribute represents the ligands bound to the site(s) of enzyme.
[8]:
# Print automatically generated names
for enzyme_form in HEX1.enzyme_module_forms:
print("Bound to sites of {0}:\n{1}\n".format(
enzyme_form.id, {
ligand.id: coeff
for ligand, coeff in enzyme_form.bound_metabolites.items()}))
Bound to sites of hex1_c:
{}
Bound to sites of hex1_A_c:
{'atp_c': 1}
Bound to sites of hex1_G_c:
{'glc__D_c': 1}
Bound to sites of hex1_AG_c:
{'glc__D_c': 1, 'atp_c': 1}
Bound to sites of hex1_G_CI_c:
{'glc__D_c': 1, '_23dpg_c': 1}
Bound to sites of hex1_A_PI_c:
{'adp_c': 1}
Bound to sites of hex1_G_PI_c:
{'g6p_c': 1}
Setting the bound_metabolites
attribute upon creation allow the formula
and charge
attributes of the various forms also to be set while ensuring mass and charge balancing is maintained. Note that the enzyme is represented as a moiety, and the ligands bound to the enzyme are represented in the chemical formula.
[9]:
# Get the elemental matrix for the enzyme
df = HEX1.get_elemental_matrix(array_type="DataFrame")
# Use iloc to only look at EnzymeModuleForms
df.iloc[:, 6:]
[9]:
hex1_c | hex1_A_c | hex1_G_c | hex1_AG_c | hex1_G_CI_c | hex1_A_PI_c | hex1_G_PI_c | |
---|---|---|---|---|---|---|---|
C | 0.0 | 10.0 | 6.0 | 16.0 | 9.0 | 10.0 | 6.0 |
H | 0.0 | 12.0 | 12.0 | 24.0 | 15.0 | 12.0 | 11.0 |
O | 0.0 | 13.0 | 6.0 | 19.0 | 16.0 | 10.0 | 9.0 |
P | 0.0 | 3.0 | 0.0 | 3.0 | 2.0 | 2.0 | 1.0 |
N | 0.0 | 5.0 | 0.0 | 5.0 | 0.0 | 5.0 | 0.0 |
S | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
q | 0.0 | -4.0 | 0.0 | -4.0 | -5.0 | -3.0 | -2.0 |
[HEX] | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 |
Setting the name
argument as “automatic” in the EnzymeModule.make_enzyme_module_form()
method causes a name for the EnzymeModuleForm
to be generated based on the metabolites in the bound_metabolites
attribute.
[10]:
# Print automatically generated names
for enzyme_form in HEX1.enzyme_module_forms:
print(enzyme_form.name)
HEX1
HEX1-atp complex
HEX1-glc__D complex
HEX1-glc__D-atp complex
HEX1-glc__D-_23dpg complex
HEX1-adp complex
HEX1-g6p complex
The categories
argument allows for EnzymeModuleForm
objects to be placed into cobra.Group
objects representing those categories. As with the ligands, the categorized enzyme module forms are returned in a DictList
of Group
objects by the enzyme_module_forms_categorized
attribute.
[11]:
for group in HEX1.enzyme_module_forms_categorized:
print("{0}: {1}".format(
group.id, str(sorted([m.id for m in group.members]))))
Active: ['hex1_AG_c', 'hex1_A_c', 'hex1_G_c', 'hex1_c']
Inhibited: ['hex1_A_PI_c', 'hex1_G_CI_c', 'hex1_G_PI_c']
Alternatively, the enzyme_module_forms_categorized
attribute can be set using a dict
:
[12]:
HEX1.enzyme_module_forms_categorized = {
"competitively_inhibited": hex1_G_CI_c}
for group in HEX1.enzyme_module_forms_categorized:
print("{0}: {1}".format(
group.id, str(sorted([m.id for m in group.members]))))
Active: ['hex1_AG_c', 'hex1_A_c', 'hex1_G_c', 'hex1_c']
Inhibited: ['hex1_A_PI_c', 'hex1_G_CI_c', 'hex1_G_PI_c']
competitively_inhibited: ['hex1_G_CI_c']
6.1.3. Defining enzyme module reactions
The next step is to define all of the reaction steps that represent the catalytic mechanism and regulation of the enzyme module. These reactions are represented as EnzymeModuleReaction
objects.
The EnzymeModuleReaction
object inherits from the MassReaction
and is treated like any other reaction in the model. Like the make_enzyme_module_form()
method, the make_enzyme_module_reaction()
method allows for the creation of an EnzymeModuleReaction
object while assigning categories for the EnzymeModuleReaction
in the process.
Species that exist in the model can also be added to the reaction by providing a dictionary of metabolites and their stoichiometric coefficients to the metabolites_to_add
argument.
[13]:
HEX1_1 = HEX1.make_enzyme_module_reaction(
"HEX1_1",
name="Automatic",
subsystem="Glycolysis",
reversible=True,
categories="product_inhibition",
metabolites_to_add={
"hex1_c": -1,
"adp_c": -1,
"hex1_A_PI_c": 1})
HEX1_2 = HEX1.make_enzyme_module_reaction(
"HEX1_2",
name="Automatic",
subsystem="Glycolysis",
reversible=True,
categories="product_inhibition",
metabolites_to_add={
"hex1_c": -1,
"g6p_c": -1,
"hex1_G_PI_c": 1})
HEX1_3 = HEX1.make_enzyme_module_reaction(
"HEX1_3",
name="Automatic",
subsystem="Glycolysis",
reversible=True,
categories="glc__D_c_binding",
metabolites_to_add={
"hex1_c": -1,
"glc__D_c": -1,
"hex1_G_c": 1})
HEX1_4 = HEX1.make_enzyme_module_reaction(
"HEX1_4",
name="Automatic",
subsystem="Glycolysis",
reversible=True,
categories="atp_c_binding",
metabolites_to_add={
"hex1_c": -1,
"atp_c": -1,
"hex1_A_c": 1})
HEX1_5 = HEX1.make_enzyme_module_reaction(
"HEX1_5",
name="Automatic",
subsystem="Glycolysis",
reversible=True,
categories="competitive_inhibition",
metabolites_to_add={
"hex1_G_c": -1,
"_23dpg_c": -1,
"hex1_G_CI_c": 1})
HEX1_6 = HEX1.make_enzyme_module_reaction(
"HEX1_6",
name="Automatic",
subsystem="Glycolysis",
reversible=True,
categories="atp_c_binding",
metabolites_to_add={
"hex1_G_c": -1,
"atp_c": -1,
"hex1_AG_c": 1})
HEX1_7 = HEX1.make_enzyme_module_reaction(
"HEX1_7",
name="Automatic",
subsystem="Glycolysis",
reversible=True,
categories="glc__D_c_binding",
metabolites_to_add={
"hex1_A_c": -1,
"glc__D_c": -1,
"hex1_AG_c": 1})
HEX1_8 = HEX1.make_enzyme_module_reaction(
"HEX1_8",
name="Automatic",
subsystem="Glycolysis",
reversible=True,
categories="catalyzation",
metabolites_to_add={
"hex1_AG_c": -1,
"hex1_c": 1,
"adp_c": 1,
"g6p_c": 1,
"h_c": 1})
for reaction in HEX1.enzyme_module_reactions:
print(reaction)
HEX1_1: adp_c + hex1_c <=> hex1_A_PI_c
HEX1_2: g6p_c + hex1_c <=> hex1_G_PI_c
HEX1_3: glc__D_c + hex1_c <=> hex1_G_c
HEX1_4: atp_c + hex1_c <=> hex1_A_c
HEX1_5: _23dpg_c + hex1_G_c <=> hex1_G_CI_c
HEX1_6: atp_c + hex1_G_c <=> hex1_AG_c
HEX1_7: glc__D_c + hex1_A_c <=> hex1_AG_c
HEX1_8: hex1_AG_c <=> adp_c + g6p_c + h_c + hex1_c
The categories
argument allows for EnzymeModuleReactions
objects to be placed into cobra.Group
objects representing those categories. As with the ligands and enzyme forms, a DictList
of the relevant groups are returned with the enzyme_module_reactions_categorized
attribute.
[14]:
HEX1.enzyme_module_reactions_categorized
[14]:
[<Group product_inhibition at 0x7fb8b0a5dac0>,
<Group glc__D_c_binding at 0x7fb911ce0340>,
<Group atp_c_binding at 0x7fb8b0a5a850>,
<Group competitive_inhibition at 0x7fb8b0a5ad30>,
<Group catalyzation at 0x7fb8b0a47b20>]
6.1.3.1. Unifying rate parameters
For this EnzymeModule
, the reactions representing glucose binding to the enzyme and ATP binding to the enzyme have the same forward rate and equilibrium constants. Instead of defining the parameter values for each individual reaction, the unify_rate_parameters()
method can be used to create custom rate laws for the given reactions that all depend on the same rate parameters.
The unify_rate_parameters()
method takes a list of reactions and an identifier to use for the unified parameter. The enzyme_prefix
flag can be set to True
to prefix the new parameter identifier with the identifier of the EnzymeModule
, ensuring that any existing custom parameters are not overwritten.
[15]:
for ligand, pid in zip([glc__D_c, atp_c],["G", "A"]):
# Get the group of reactions corresponding to the ligand
category = "_".join((ligand.id, "binding"))
group = HEX1.enzyme_module_reactions_categorized.get_by_id(category)
# Unify the parameters
HEX1.unify_rate_parameters(
group.members, new_parameter_id=pid, enzyme_prefix=True)
# Print the new reaction rates
print("\n" + category + "\n" + "-" * len(category))
for reaction in sorted(group.members, key=attrgetter("id")):
print(reaction.id + ": " + str(reaction.rate))
glc__D_c_binding
----------------
HEX1_3: kf_HEX1_G*(glc__D_c(t)*hex1_c(t) - hex1_G_c(t)/Keq_HEX1_G)
HEX1_7: kf_HEX1_G*(glc__D_c(t)*hex1_A_c(t) - hex1_AG_c(t)/Keq_HEX1_G)
atp_c_binding
-------------
HEX1_4: kf_HEX1_A*(atp_c(t)*hex1_c(t) - hex1_A_c(t)/Keq_HEX1_A)
HEX1_6: kf_HEX1_A*(atp_c(t)*hex1_G_c(t) - hex1_AG_c(t)/Keq_HEX1_A)
6.2. Determining Enzyme Form Concentrations and Rate Constants
The next step is to solve for the steady state concentrations for the various forms of the enzyme symbolically using SymPy. Because the numerical values for the dissociation constants have been defined, these equations are solved in terms of the rate constants. The rate constants can be approximated using the total enzyme concentration as a constraint and substituted back into the equations to calculate the numerical values of the steady state concentrations.
[16]:
from sympy import Eq, Symbol, lambdify, simplify, solveset
from mass import strip_time
from mass.util.matrix import matrix_rank
6.2.1. Solving steady state concentrations symbolically
To get the symbolic solutions for the individual enzyme forms, the ODEs are first collected in a dict
. Keys are the enzyme forms, and values are their ODEs with the time dependency stripped via the strip_time
function.
[17]:
ode_dict = {
enzyme_form.id: Eq(strip_time(enzyme_form.ode), 0)
for enzyme_form in HEX1.enzyme_module_forms}
# Matrix rank of enzyme stoichiometric matrix without substrates
rank = matrix_rank(HEX1.S[6:])
print("Rank Deficiency: {0}".format(len(ode_dict) - rank))
Rank Deficiency: 1
Because the stoichiometric matrix (without ligands) has a rank deficiency of one, there is a dependent variable in the system unless another equation is added. Therefore, the completely free enzyme form is treated as the dependent variable, and all of the enzyme forms are solved in terms of the free enzyme form.
[18]:
enzyme_solutions = {}
for enzyme_form in HEX1.enzyme_module_forms:
# Skip dependent variable
if enzyme_form.id == "hex1_c":
continue
# Get the ODE for the enzyme form from the ODE dict
equation = ode_dict[enzyme_form.id]
# Solve the equation for the enzyme form, substituting
# previously found enzyme form solutions into the equation
solution = solveset(equation.subs(enzyme_solutions),
enzyme_form.id)
# Store the solution
enzyme_solutions[enzyme_form.id] = list(solution)[0]
# Substitute the new solution into existing solutions
enzyme_solutions.update({
enzyme_form: sol.subs(enzyme_solutions)
for enzyme_form, sol in enzyme_solutions.items()})
args = set()
for solution in enzyme_solutions.values():
args.update(solution.atoms(Symbol))
6.2.1.1. Defining the Rate Equation
To make up for the rank deficiency, an additional equation is needed. Typically, the rate of the enzyme is the summation of the rates for the catalyzation reaction step(s) of the enzyme. The make_enzyme_rate_equation()
method can be used to create the rate equation from a list of reactions. If use_rates=True
, the rate expressions of the reactions are added together. If update_enzyme=True
, the rate equation is set as a symbolic expression for the enzyme_rate_equation
attribute.
[19]:
# Get the catalyzation reactions
catalyzation_group = HEX1.enzyme_module_reactions_categorized.get_by_id(
"catalyzation")
HEX1.make_enzyme_rate_equation(catalyzation_group.members,
use_rates=True,
update_enzyme=True)
print(HEX1.enzyme_rate_equation)
kf_HEX1_8*(Keq_HEX1_8*hex1_AG_c(t) - adp_c(t)*g6p_c(t)*hex1_c(t))/Keq_HEX1_8
With the rate equation defined, the enzyme_rate_error()
method is used to get the equation as the difference between the flux value and the rate equation.
[20]:
enzyme_rate_equation = strip_time(HEX1.enzyme_rate_error(use_values=False))
print(enzyme_rate_equation)
v_HEX1 - kf_HEX1_8*(Keq_HEX1_8*hex1_AG_c - adp_c*g6p_c*hex1_c)/Keq_HEX1_8
The solutions for the enzyme forms are substituted into the rate equation, and the equation is solved for the free enzyme form. The solutions are subsequently updated, resulting in symbolic equations that do not depend on any enzyme form.
[21]:
# Solve for last unknown concentration symbolically
solution = solveset(enzyme_rate_equation.subs(enzyme_solutions),
"hex1_c")
# Update solution dictionary with the new solution
enzyme_solutions["hex1_c"] = list(solution)[0]
# Update solutions with free variable solutions
enzyme_solutions = {
enzyme_form: simplify(solution.subs(enzyme_solutions))
for enzyme_form, solution in enzyme_solutions.items()}
args = set()
for solution in enzyme_solutions.values():
args.update(solution.atoms(Symbol))
print(args)
{atp_c, Keq_HEX1_8, glc__D_c, kf_HEX1_G, adp_c, Keq_HEX1_5, v_HEX1, Keq_HEX1_2, kf_HEX1_A, kf_HEX1_8, Keq_HEX1_1, _23dpg_c, Keq_HEX1_G, Keq_HEX1_A, g6p_c}
Numerical values for known quantities are substituted into the equations. For this EnzymeModule
of Hexokinase, the following dissociation constants are used:
A value of \(K_{\text{HEX1}}= 313.12\) is used for the catalyzation step. Note that the inverse of the dissociation constant is used for reactions that form complexes.
[22]:
numerical_values = {
"Keq_HEX1_1": 1,
"Keq_HEX1_2": 1 / 66.67,
"Keq_HEX1_G": 1 / 0.038,
"Keq_HEX1_A": 1 / 2.06,
"Keq_HEX1_5": 1 / 5.5,
"Keq_HEX1_8": 313.12}
# Update the model with the parameters
HEX1.update_parameters(numerical_values)
The ligand concentrations and the rate for the enzyme are extracted from the merged glycolysis and hemoglobin model.
[23]:
# Get steady state flux for EnzymeModule
HEX1.enzyme_rate = glyc_hb.reactions.get_by_id("HEX1").steady_state_flux
numerical_values[HEX1.enzyme_flux_symbol_str] = HEX1.enzyme_rate
# Get the ligand concentrations
for met in HEX1.enzyme_module_ligands:
concentration = glyc_hb.metabolites.get_by_id(met.id).initial_condition
# Set the ligand initial condition and add to numercal values dictionary
met.initial_condition = concentration
numerical_values[met.id] = concentration
The numerical values are substituted into the symbolic equations, resulting in the steady state concentrations that depend only on the rate constants.
[24]:
enzyme_solutions = {
enzyme_form: simplify(sol.subs(numerical_values))
for enzyme_form, sol in enzyme_solutions.items()}
args = set()
for solution in enzyme_solutions.values():
args.update(solution.atoms(Symbol))
print(args)
{kf_HEX1_8, kf_HEX1_G, kf_HEX1_A}
6.2.2. Approximating Rate Constants
To determine the set of rate constants for the enzyme module, the absolute error between the total hexokinase concentration value (found in literature) and the computed hexokinase concentration is minimized. For this example, the minimize()
function of the SciPy package is utilized to find a feasible set of rate constants.
[25]:
from scipy.optimize import minimize
The objective function for the minimization is first made symbolically. The enzyme_total_symbol_str
property can be used to represent the total enzyme concentration, while the enzyme_concentration_total_equation
property creates a symbolic expression for the sum of all enzyme forms.
[26]:
enzyme_total_error = abs(
Symbol(HEX1.enzyme_total_symbol_str)
- strip_time(HEX1.enzyme_concentration_total_equation))
print(enzyme_total_error)
Abs(-HEX1_Total + hex1_AG_c + hex1_A_PI_c + hex1_A_c + hex1_G_CI_c + hex1_G_PI_c + hex1_G_c + hex1_c)
The enzyme_concentration_total
attribute stores the total amount of enzyme in the model and substituted into the expression. The total HEX1 concentration is \(24 * 10^{-6} \text{mM}\).
[27]:
HEX1.enzyme_concentration_total = 24e-6
enzyme_total_error = enzyme_total_error.subs({
HEX1.enzyme_total_symbol_str: HEX1.enzyme_concentration_total})
print(enzyme_total_error)
Abs(hex1_AG_c + hex1_A_PI_c + hex1_A_c + hex1_G_CI_c + hex1_G_PI_c + hex1_G_c + hex1_c - 2.4e-5)
Finally, the symbolic equations for the enzyme forms are substituted into the enzyme total error equation, resulting in an expression that represents the objective function with the only unknown variables being rate constants. The lambdify()
function of the SymPy package converts the symbolic objective into a lambda function that can be used with the minimize()
function of SciPy.
[28]:
enzyme_total_error = simplify(enzyme_total_error.subs(enzyme_solutions))
# Sort the arguments to ensure input format remains consistent
args = sorted(list(map(str, args)))
# Use lambdify to make objective function as a lambda function
obj_fun = lambda x: lambdify(args, enzyme_total_error)(*x)
The minimize()
function is now used to approximate the rate constants. The optimization problems for enzyme rate constants are typically nonlinear, and require nonlinear optimization routines to find feasible solutions.
[29]:
# Minimize the objective function, initial guess based on publication values
initial_guess = [1e8, 9376585, 52001]
variable_bounds = ((0, 1e9), (0, 1e9), (0, 1e9))
solution = minimize(obj_fun, x0=initial_guess,
method="trust-constr",
bounds=variable_bounds)
# Map solution array to variables
rate_constants = dict(zip(args, solution.x))
print(rate_constants)
{'kf_HEX1_8': 100000000.0025878, 'kf_HEX1_A': 9376585.030755484, 'kf_HEX1_G': 52006.59981223971}
Because the rate constants associated with the inhibition of the enzyme forms are not necessary for computing the concentrations, a rapid binding assumption is made for the inhibition reactions. Therefore, a large number is set for the rate constants. The parameters are set using the update_parameters()
method.
[30]:
rate_constants["kf_HEX1_1"] = 1e6
rate_constants["kf_HEX1_2"] = 1e6
rate_constants["kf_HEX1_5"] = 1e6
HEX1.update_parameters(rate_constants)
6.2.3. Calculating numerical values for concentrations
Once the rate constants have been estimated, they are substituted back into the symbolic concentration equations in order to obtain their numerical values.
[31]:
for enzyme_form, solution in enzyme_solutions.items():
# Get the enzyme form object, determine the steady state concentration
enzyme_form = HEX1.enzyme_module_forms.get_by_id(enzyme_form)
enzyme_form.initial_condition = float(solution.subs(rate_constants))
print("{0}: {1:e}".format(enzyme_form.id,
enzyme_form.initial_condition))
hex1_A_c: 9.401421e-06
hex1_G_c: 5.718872e-08
hex1_AG_c: 1.174630e-08
hex1_G_CI_c: 3.223364e-08
hex1_A_PI_c: 3.519706e-06
hex1_G_PI_c: 8.847367e-09
hex1_c: 1.213692e-05
6.2.3.1. Error values
As a quality assurance check, the enzyme_concentration_total_error()
method can be used to get the error between the enzyme_concentration_total
attribute and the sum of the enzyme form concentrations. A positive value indicates the enzyme_concentration_total
attribute is greater than the sum of the individual enzyme form concentrations that were computed.
[32]:
print("Total Enzyme Concentration Error: {0}".format(
HEX1.enzyme_concentration_total_error(use_values=True)))
Total Enzyme Concentration Error: -1.1680622689093244e-06
Similarly, the error between the enzyme_rate
attribute and the computed value from the enzyme_rate_equation
can be also checked using the enzyme_rate_error()
method, in which a positive value indicates that the enzyme_rate
attribute is greater than the value computed when using the rate equation.
[33]:
print("Enzyme Rate Error: {0}".format(
HEX1.enzyme_rate_error(use_values=True)))
Enzyme Rate Error: 4.440892098500626e-16
6.3. Adding EnzymeModules to Models
With the EnzymeModule
built, it can be integrated into a larger network and simulated. To add an EnzymeModule
to an existing MassModel
, the merge()
method is used. After merging, the remove_reactions()
method is used to remove the reaction replaced with the enzyme module. The EnzymeModule
should always be merged into the MassModel
as demonstrated below:
[34]:
glyc_hb_HEX1 = glyc_hb.merge(HEX1, inplace=False)
glyc_hb_HEX1.remove_reactions([
glyc_hb_HEX1.reactions.get_by_id("HEX1")])
All objects, numerical values, and certain attributes of the EnzymeModule
are transferred into the MassModel
upon merging. This includes all enzyme forms, reactions steps, initial conditions, rate parameters, and category groups.
[35]:
glyc_hb_HEX1
[35]:
Name | Glycolysis_Hemoglobin_HEX1 |
Memory address | 0x07fb911fdf2b0 |
Stoichiometric Matrix | 35x37 |
Matrix Rank | 32 |
Number of metabolites | 35 |
Initial conditions defined | 35/35 |
Number of reactions | 37 |
Number of genes | 0 |
Number of enzyme modules | 1 |
Number of groups | 12 |
Objective expression | 0 |
Compartments | Cytosol |
6.3.1. The EnzymeModuleDict object
During the merge process, an EnzymeModuleDict
is created from the EnzymeModule
and added to the MassModel.enzyme_modules
attribute.
[36]:
print(glyc_hb_HEX1.enzyme_modules)
HEX1_dict = glyc_hb_HEX1.enzyme_modules.get_by_id("HEX1")
HEX1_dict
[<EnzymeModuleDict HEX1 at 0x7fb911e56940>]
[36]:
Name | HEX1 |
Memory address | 0x07fb911e56940 |
Stoichiometric Matrix | 13x8 |
Matrix Rank | 7 |
Subsystem | Glycolysis |
Number of Ligands | 6 |
Number of EnzymeForms | 7 |
Number of EnzymeModuleReactions | 8 |
Enzyme Concentration Total | 2.4e-05 |
Enzyme Net Flux | 1.12 |
The EnzymeModuleDict
inherits from an OrderedDict
, thereby inheriting ordered dictionary methods such as keys()
:
[37]:
print("\n".join(HEX1_dict.keys()))
id
name
subsystem
enzyme_module_ligands
enzyme_module_forms
enzyme_module_reactions
enzyme_module_ligands_categorized
enzyme_module_forms_categorized
enzyme_module_reactions_categorized
enzyme_concentration_total
enzyme_rate
enzyme_concentration_total_equation
enzyme_rate_equation
S
model
The EnzymeModuleDict
stores several of the enzyme-specific attributes so that they are still accessible after integrating the enzyme module into a larger network. The keys of the EnzymeModuleDict
also can be treated as attribute accessors:
[38]:
print("Enzyme Rate:\n{0} = {1}".format(
HEX1_dict["enzyme_rate"], # Returned using dict key
HEX1_dict.enzyme_rate_equation # Returned using attribute accessor
))
Enzyme Rate:
1.12 = kf_HEX1_8*(Keq_HEX1_8*hex1_AG_c(t) - adp_c(t)*g6p_c(t)*hex1_c(t))/Keq_HEX1_8
6.3.2. Steady State Validation
The last step is to ensure that a steady state is reached with the completed enzyme module within a larger network context.
[39]:
import matplotlib.pyplot as plt
from mass import Simulation
from mass.visualization import plot_time_profile
Here, the model is simulated, and the enzyme’s ability to reach steady state is graphically verified:
[40]:
# Setup simulation object
sim = Simulation(glyc_hb_HEX1, verbose=True)
# Simulate from 0 to 1000 with 10001 points in the output
conc_sol, flux_sol = sim.simulate(
glyc_hb_HEX1, time=(0, 1e3))
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 4))
plot_time_profile(
conc_sol, observable=HEX1_dict.enzyme_module_forms, ax=ax,
xlim=(1e-3, 1e3),
legend="right outside", plot_function="loglog",
xlabel="Time [hr]", ylabel="Concentration [mM]",
title="TIme profile of Concentrations for Enzyme Forms");
Successfully loaded MassModel 'Glycolysis_Hemoglobin_HEX1' into RoadRunner.
The plot shows that the enzyme can reach a steady state when integrated into a larger network, meaning the enzyme module that represents hexokinase in this system is complete!
6.4. Additional Examples
For additional examples of analyzing and visualizing systems with enzyme modules, see the following:
\(^1\) Procedure outlined in [DZK+16]
\(^2\) Hexokinase based on [YAHP18], [DZK+16], and [MBK99]
\(^3\) Glycolysis model based on [YAHP18] and Chapter 10 of [Pal11]