!pip install masspy

import numpy as np import pandas as pd from scipy import optimize from cobra import DictList from mass import ( MassConfiguration, MassMetabolite, MassModel, MassReaction, MassModel, Simulation, MassSolution, strip_time, Unit, UnitDefinition, qcqa_model) from mass.util.matrix import left_nullspace, nullspace, matrix_rank from mass.util.expressions import Keq2k, k2Keq, strip_time from mass.util.qcqa import qcqa_model from mass.enzyme_modules import EnzymeModule from mass.core.units import print_defined_unit_values from operator import attrgetter from os import path import matplotlib as mpl import matplotlib.pyplot as plt import sympy as sym from mass.io import json, sbml from mass.visualization import ( plot_time_profile, plot_phase_portrait, plot_tiled_phase_portraits)

# Configuration setup mass_config = MassConfiguration() mass_config.irreversible_Keq = float("inf")

# Visualization setup pd.set_option("display.max_rows", 100) pd.set_option("display.max_columns", 100) pd.set_option('display.max_colwidth', None) pd.options.display.float_format = '{:,.3f}'.format S_FONT = {"size": "small"} L_FONT = {"size": "large"} INF = float("inf")

model = MassModel("Urea_Cycle")

""" Each metabolite in the network is created with an ID, a name, its formula, its charge, and the compartment it is in (mitochrondia, cytosol or both). """ hco3_m = MassMetabolite( "hco3_m", name="Carbonic acid", formula="HCO3", charge=-1, compartment="m") nh4_m = MassMetabolite( "nh4_m", name="Ammonium", formula="NH4", charge=1, compartment="m") h = MassMetabolite( "h", name="Hydrogen", formula="H", charge=1, compartment="both") atp = MassMetabolite( "atp", name="ATP", formula="C10H12N5O13P3", charge=-4, compartment="both") adp_m = MassMetabolite( "adp_m", name="ADP", formula="C10H12N5O10P2", charge=-3, compartment="m") pi_m = MassMetabolite( "pi_m", name="Phosphate", formula="PO4", charge=-3, compartment="m") carb_phos_m = MassMetabolite( "carb_phos_m", name="Carbamoyl phosphate", formula="NH2CO2PO3", charge=-2, compartment="m") ornithine = MassMetabolite( "ornithine", name="Ornithine", formula="C5H13N2O2", charge=1, compartment="both") citrulline = MassMetabolite( "citrulline", name="Citrulline", formula="C6H13N3O3", charge=0, compartment="both") ppi_c = MassMetabolite( "ppi_c", name="Diphosphate", formula="P2O7", charge=-4, compartment="c") amp_c = MassMetabolite( "amp_c", name="AMP", formula="C10H12N5O7P", charge=-2, compartment="c") aspartate_c = MassMetabolite( "aspartate_c", name="Aspartate", formula="C4H6NO4", charge=-1, compartment="c") argi_suc_c = MassMetabolite( "argi_suc_c", name="Argigino succinate", formula="C10H17N4O6", charge=-1, compartment="c") fumarate_c = MassMetabolite( "fumarate_c", name="Fumarate", formula="C4H2O4", charge=-2, compartment="c") arginine_c = MassMetabolite( "arginine_c", name="Arginine", formula="C6H15N4O2", charge=1, compartment="c") h2o_c = MassMetabolite( "h2o_c", name="Water", formula="H2O", charge=0, compartment="c") urea_c = MassMetabolite( "urea_c", name="Urea", formula="CH4N2O", charge=0, compartment="c")

# Define main reactions in the network enter = MassReaction("CPS", name="Entering", reversible=False) condense_1 = MassReaction("OTC", name="Condensation 1", reversible=False) condense_2 = MassReaction("ASD", name="Condensation 2", reversible=False) cleave_1 = MassReaction("ALD", name="Cleavage 1", reversible=False) cleave_2 = MassReaction("AG", name="Cleavage 2", reversible=False) # Add metabolites # the value of each metabolite shows whether it is a reactant (negative) or a product (positive), # and its coefficient enter.add_metabolites({ nh4_m: -1, hco3_m: -1, atp: -2, h: 3, pi_m: 1, adp_m: 2, carb_phos_m: 1 }) condense_1.add_metabolites({ carb_phos_m: -1, ornithine: -1, pi_m: 1, h: 2, citrulline: 1 }) condense_2.add_metabolites({ citrulline: -1, aspartate_c: -1, atp: -1, h: 2, ppi_c: 1, amp_c: 1, argi_suc_c: 1 }) cleave_1.add_metabolites({ argi_suc_c: -1, fumarate_c: 1, arginine_c: 1 }) cleave_2.add_metabolites({ arginine_c: -1, h2o_c: -1, ornithine: 1, urea_c: 1 }) # Add the reactions to the model model.add_reactions([enter, condense_1, condense_2, cleave_1, cleave_2])

# Define and add boundary reactions # for example, there is a source for hco3 hco3_in = model.add_boundary( model.metabolites.hco3_m, boundary_type="demand", reaction_id="hco3_in") # by default, it is consumed, so we have to reverse it hco3_in.reverse_stoichiometry(inplace=True) nh4_in = model.add_boundary( model.metabolites.nh4_m, boundary_type="demand", reaction_id="nh4_in") nh4_in.reverse_stoichiometry(inplace=True) atp_in = model.add_boundary( model.metabolites.atp, boundary_type="demand", reaction_id="atp_in") atp_in.reverse_stoichiometry(inplace=True) aspartate_in = model.add_boundary( model.metabolites.aspartate_c, boundary_type="demand", reaction_id="aspartate_in") aspartate_in.reverse_stoichiometry(inplace=True) h2o_in = model.add_boundary( model.metabolites.h2o_c, boundary_type="demand", reaction_id="h2o_in") h2o_in.reverse_stoichiometry(inplace=True) h_out = model.add_boundary( model.metabolites.h, boundary_type="demand", reaction_id="h_out") adp_out = model.add_boundary( model.metabolites.adp_m, boundary_type="demand", reaction_id="adp_out") pi_out = model.add_boundary( model.metabolites.pi_m, boundary_type="demand", reaction_id="pi_out") amp_out = model.add_boundary( model.metabolites.amp_c, boundary_type="demand", reaction_id="amp_out") ppi_out = model.add_boundary( model.metabolites.ppi_c, boundary_type="demand", reaction_id="ppi_out") fumarate_out = model.add_boundary( model.metabolites.fumarate_c, boundary_type="demand", reaction_id="fumarate_out") urea_out = model.add_boundary( model.metabolites.urea_c, boundary_type="demand", reaction_id="urea_out")

# Check the number of reactions and metabolites print(f"Number of metabolites: {len(model.metabolites)}") print(f"Number of all reactions: {len(model.reactions)}") print(f"Number of boundary reactions: {len(model.boundary)}")

# Check all the boundary reactions print('All the boundary reactions') for boundary in model.boundary: print(boundary) # Check all the reactions print('All the main reactions') for reaction in model.reactions[:5]: print(reaction)

np.random.seed(123)

# Add the initial concentrations of all metabolites for metabolite in model.metabolites: metabolite.ic = np.random.uniform(0, 10) h2o_c.ic = 55000 hco3_m.ic = 24.5 nh4_m.ic = 0.0385 h.ic = 10 atp.ic = 1.077 adp_m.ic = 0.16 pi_m.ic = 1.137 carb_phos_m.ic = 0.029 ornithine.ic = 0.013 citrulline.ic = 0.02875 ppi_c.ic = 0.0018 amp_c.ic = 0.019 aspartate_c.ic = 0.0155 argi_suc_c.ic = 0.00175 fumarate_c.ic = 0.000945 arginine_c.ic = 0.07683 urea_c.ic = 1.682 # Show all the metabolites and their initial concentrations metabolite_ids = [m.id for m in model.metabolites] table_1 = pd.DataFrame( np.array([metabolite_ids, [met.name for met in model.metabolites], [model.initial_conditions[met] for met in model.metabolites]]).T, index=[i for i in range(1, len(metabolite_ids) + 1)], columns=["Abbreviations", "Species", "Initial Concentration"]) table_1

# hco3_in.kf = 2.6 # nh4_in.kf = 6 # atp_in.kf = 1.2 # aspartate_in.kf = 3.5 # h2o_in.kf = 3.2 # h_out.kf = 4.3 # adp_out.kf = 5.2 # pi_out.kf = 3.1 # amp_out.kf = 2.7 # ppi_out.kf = 4.1 # fumarate_out.kf = 1.2 # urea_out.kf = 2.8 for reaction in model.reactions: reaction.kf = np.random.uniform(200, 2000) enter.forward_rate_constant = 1.213*1000 condense_1.kf = 2.3*1000 condense_2.kf = 0.1295*1000 cleave_1.kf = 10.3*1000 cleave_2.kf = 0.034*1000 for reaction in model.reactions: print(reaction.parameters)

# Show the stoichiometric matrix model.update_S(array_type="DataFrame", dtype=np.int_, update_model=True) table_2 = model.S table_2

print('Dimensions of matrix: ', model.S.shape) print('Rank of matrix: ', matrix_rank(model.S))

# Show the elemental matrix table_3 = model.get_elemental_matrix( array_type="DataFrame", dtype=np.int64) table_3

# Check if the model is elementally and charge balanced table_4 = model.get_elemental_charge_balancing( array_type="DataFrame", dtype=np.int64) table_4

# Get the null space ns = nullspace(model.S) # Divide by the minimum and round to nearest integer ns = np.rint(ns / np.min(ns[np.nonzero(ns)])) # Show the nullspace table_5 = pd.DataFrame( ns, index=model.reactions.list_attr("id"), columns=["Pathway 1"], dtype=np.int_).T table_5

# Check if the found pathway is balanced table_6 = model.get_elemental_charge_balancing( array_type="DataFrame").dot(ns).astype(np.int64) table_6.columns = ['Pathway 1'] table_6

# Get the left nullspace lns = left_nullspace(model.S, rtol=1e-10) # Iterate through left nullspace, dividing by the smallest value in each row. for i, row in enumerate(lns): minval = np.min(abs(row[np.nonzero(row)])) new_row = np.array(row/minval) # Round to ensure the left nullspace is composed of only integers lns[i] = np.array([round(value) for value in new_row]) # Ensure positive stoichiometric coefficients if all are negative for i, space in enumerate(lns): lns[i] = np.negative(space) if all([num <= 0 for num in space]) else space # Show the left nullspace lns = lns.astype(np.int64) table_7 = pd.DataFrame(lns, columns=metabolite_ids) table_7

model.compartments = {"c": "cytosol", "m": "mitochondria", "both": "both"} model.compartments

# Check the default units print_defined_unit_values("BaseUnitKinds")

# Create concentration unit millimole = Unit(kind="mole", exponent=1, scale=-3, multiplier=1) per_liter = Unit(kind="liter", exponent=-1, scale=1, multiplier=1) concentration_unit = UnitDefinition(id="mM", name="Millimolar", list_of_units=[millimole, per_liter]) print('Concentration unit') print("{0}:\n{1!r}\n{2!r}".format( concentration_unit.name, *concentration_unit.list_of_units)) # Create time unit time_unit = UnitDefinition(id="min", name="Minute") time_unit.create_unit(kind="second", exponent=1, scale=1, multiplier=60) print('Time unit') print(time_unit) print(time_unit.list_of_units) # Add units to the model model.add_units([concentration_unit, time_unit]) model.units

model

# Check the rates of reactions for reaction, rate in model.rates.items(): print(f"{reaction.id}: {rate}")

# Check the ODE of each metabolite for metabolite, ode in model.odes.items(): print(f"{metabolite.id}: {ode}")

qcqa_model(model, parameters=True, concentrations=True)

# Set independent fluxes to determine steady state flux vector independent_fluxes = { model.reactions.nh4_in: 2.1} # Compute steady state fluxes ssfluxes = model.compute_steady_state_fluxes( ns.T, independent_fluxes, update_reactions=True) # Show the steady ste fluxes reaction_ids = [r.id for r in model.reactions] table_8 = pd.DataFrame(list(ssfluxes.values()), index=reaction_ids, columns=["v_stst"]).T table_8

# Visualize steady state fluxes # Set up the figure fig_1, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 5)) # Define indicies for bar chart indicies = np.arange(len(reaction_ids))+0.5 # Define colors to use c = plt.cm.coolwarm(np.linspace(0, 1, len(reaction_ids))) # Plot bar chart ax.bar(indicies, ssfluxes.values(), width=0.8, color=c); ax.set_xlim([0, len(reaction_ids)]) # Set labels and adjust ticks ax.set_xticks(indicies) ax.set_xticklabels(reaction_ids, rotation="vertical") ax.set_ylabel("Fluxes (mM/hr)", L_FONT) ax.set_title("Steady State Fluxes", L_FONT) fig_1.tight_layout()

# Check if this is really the steady state flux pd.DataFrame( model.S.dot(np.array(list(ssfluxes.values()))), index=metabolite_ids, columns=["S x v_stst"], dtype=np.int64).T

percs = model.calculate_PERCs(update_reactions=True) percs

for r in model.reactions: print(r.parameters)

sim0 = Simulation(reference_model=model) sim0.steady_state_solver.maximum_iterations = 5000 conc_sol, flux_sol = sim0.find_steady_state(model, strategy="simulate", update_values=True, verbose=True) model.initial_conditions # Same object as reference model in Simulation

for metabolite, solution in conc_sol.items(): print("{0}: {1}".format(metabolite, solution))

# Simulate the model t0, tf = (0, 1e3) sim = Simulation(model, verbose=True)

# Visualize the concentrations over time conc_sol_ss, flux_sol_ss = sim.simulate( model, time=(t0, tf, tf*10 + 1)) conc_sol_ss.view_time_profile()

# Plot time profile plot_time_profile( conc_sol_ss, legend="right outside", xlabel="Time", ylabel="Concentration", title=("Concentration Solutions", {"size": "large"}))

# Plot time profile for concentrations on log scale plot_time_profile( conc_sol_ss, legend="right outside", plot_function="semilogx", xlabel="Time", ylabel="Concentration", title=("Concentration Solutions", {"size": "large"}))

# Plot time profil for fluxes on log scale plot_time_profile( flux_sol_ss, legend="right outside", plot_function="semilogx", xlabel="Time [hr]", ylabel="Flux [mM/hr]", title=("Flux Profile", L_FONT));

# Plot time profile of ammonium and urea plot_time_profile( conc_sol_ss, observable=["nh4_m", "carb_phos_m", "ornithine", "citrulline", "aspartate_c", "argi_suc_c", "urea_c"], legend="right outside", plot_function="semilogx", xlabel="Time", ylabel="Concentration", title=("Concentration Solutions", {"size": "large"}))

# Plot phase portrait of ammonium and urea plot_phase_portrait( conc_sol_ss, x="nh4_m", y="urea_c", xlabel="Ammonia", ylabel="Urea", annotate_time_points=[t0, 1e-1, 1e0, 1e1, 1e3, tf], annotate_time_points_color= [ "red", "green", "purple", "yellow", "cyan", "blue"], annotate_time_points_legend="lower outside")

# Plot tiled phase portraits of all metabolites # plot_tiled_phase_portraits(conc_sol_ss, # annotate_time_points_legend="right outside")

CPS = EnzymeModule("CPS", name="Ammonia-dependent carbamoyl-phosphate synthase", subsystem="Urea_Cycle")

# Add the metabolites to the EnzymeModule CPS.add_metabolites([nh4_m, hco3_m, atp, h, pi_m, adp_m, carb_phos_m]) # Access DictList of ligands and print print("All {0} Ligands: {1}".format( CPS.id, "; ".join([m.id for m in CPS.enzyme_module_ligands]))) CPS.enzyme_module_ligands_categorized = { "Substrates": [nh4_m, hco3_m], "Cofactors": atp, "Products": [h, pi_m, adp_m, carb_phos_m]} # Access DictList of ligands and print print("All {0} ligands ({1} total):\n{2}\n".format( CPS.id, len(CPS.enzyme_module_ligands), str([m.id for m in CPS.enzyme_module_ligands]))) # Access categorized attribute for ligands and print for group in CPS.enzyme_module_ligands_categorized: print("{0}: {1}".format( group.id, str([m.id for m in group.members])))

# Make enzyme module forms per number of bound activators (Up to 4 Total) CPS.make_enzyme_module_form( "cps_m", name="Automatic", categories=("Relaxed"), compartment="m"); CPS.make_enzyme_module_form( "cps_N_m", name="Automatic", categories=("Relaxed"), bound_metabolites={nh4_m: 1}, compartment="m"); CPS.make_enzyme_module_form( "cps_NH_m", name="Automatic", categories=("Relaxed"), bound_metabolites={nh4_m: 1, hco3_m: 1}, compartment="m"); CPS.make_enzyme_module_form( "cps_NHA_m", name="Automatic", categories=("Relaxed"), bound_metabolites={nh4_m: 1, hco3_m: 1, atp: 1, atp: 2}, compartment="m"); # Access DictList of enzyme module forms and print print("All {0} enzyme module forms ({1} total):\n{2}\n".format( CPS.id, len(CPS.enzyme_module_forms), str([m.id for m in CPS.enzyme_module_forms]))) # Access categorized attribute for enzyme module forms and print for group in CPS.enzyme_module_forms_categorized: print("{0}: {1}\n".format( group.id, str(sorted([m.id for m in group.members]))))

CPS.make_enzyme_module_reaction( "CPS_1".format(i), name="Automatic", subsystem="Urea_Cycle", reversible=True, categories="nh4_m_binding", metabolites_to_add={ "cps_m": -1, "nh4_m": -1, "cps_N_m": 1}) CPS.make_enzyme_module_reaction( "CPS_2", name="Automatic", subsystem="Urea_Cycle", reversible=True, categories="hco3_m_binding", metabolites_to_add={ "cps_N_m": -1, "hco3_m": -1, "cps_NH_m": 1}) CPS.make_enzyme_module_reaction( "CPS_3", name="Automatic", subsystem="Urea_Cycle", reversible=True, categories="atp_binding", metabolites_to_add={ "cps_NH_m": -1, "atp": -2, "cps_NHA_m": 1}) CPS.make_enzyme_module_reaction( "CPS_4", name="Automatic", subsystem="Urea_Cycle", reversible=False, categories="catalyzation", metabolites_to_add={ "cps_NHA_m": -1, "cps_m": 1, "h": 2, "pi_m": 1, "adp_m": 2, "carb_phos_m": 1}) # Access DictList of enzyme module reactions and print print("All {0} enzyme module reactions ({1} total):\n{2}\n".format( CPS.id, len(CPS.enzyme_module_reactions), str([m.name for m in CPS.enzyme_module_reactions]))) # Access categorized attribute for enzyme module reactions and print for group in CPS.enzyme_module_reactions_categorized: print("{0}: {1}\n".format( group.id, str(sorted([m.id for m in group.members]))))

abbreviations = ["N", "H", "A"] ligands = [nh4_m, hco3_m, atp, pi_m, adp_m, h, carb_phos_m] for met, unified_id in zip(ligands, abbreviations): category = {"N": "binding", "H": "binding", "A": "binding" }[unified_id] group = CPS.enzyme_module_reactions_categorized.get_by_id( "_".join((met.id, category))) reactions = sorted(group.members, key=attrgetter("id")) CPS.unify_rate_parameters(reactions, unified_id, rate_type=2, enzyme_prefix=True) CPS.unify_rate_parameters( CPS.enzyme_module_reactions_categorized.get_by_id("catalyzation").members, "CPS") # Update rate laws to be in terms of kf and Keq CPS.custom_rates.update(k2Keq(CPS.custom_rates)) # Access categorized attribute for enzyme module reactions and print for group in CPS.enzyme_module_reactions_categorized: header = "Category: " + group.id print("\n" + header + "\n" + "-" * len(header)) for reaction in sorted(group.members, key=attrgetter("id")): print(reaction.id + ": " + str(reaction.rate))

# Make a dictionary of ODEs and lump ligands into rate parameters by giving them a value of 1 ode_dict = {} lump_ligands = {sym.Symbol(met.id): 1 for met in CPS.enzyme_module_ligands} for enzyme_module_form in CPS.enzyme_module_forms: symbol_key = sym.Symbol(enzyme_module_form.id) ode = sym.Eq(strip_time(enzyme_module_form.ode), 0) ode_dict[symbol_key] = ode.subs(lump_ligands) rank = matrix_rank(CPS.S[7:]) print("Rank Deficiency: {0}".format(len(ode_dict) - rank))

reactions = CPS.enzyme_module_reactions_categorized.get_by_id( "catalyzation").members CPS.make_enzyme_rate_equation( reactions, use_rates=True, update_enzyme=True) sym.pprint(CPS.enzyme_rate_equation)

# Get enzyme module forms enzyme_module_forms = CPS.enzyme_module_forms.copy() # Reverse list for increased performance (due to symmetry assumption) # by solving for the most activated/inhibitors bound first. enzyme_module_forms.reverse() enzyme_solutions = {} for enzyme_module_form in enzyme_module_forms: # Skip dependent variable if "cps_m" == str(enzyme_module_form): continue enzyme_module_form = sym.Symbol(enzyme_module_form.id) # Susbtitute in previous solutions and solve for the enzyme module form, equation = ode_dict[enzyme_module_form] sol = sym.solveset(equation.subs(enzyme_solutions), enzyme_module_form) enzyme_solutions[enzyme_module_form] = list(sol)[0] # Update the dictionary of solutions with the solutions enzyme_solutions.update({ enzyme_module_form: sol.subs(enzyme_solutions) for enzyme_module_form, sol in enzyme_solutions.items()}) args = set() for sol in enzyme_solutions.values(): args.update(sol.atoms(sym.Symbol)) print(args)

enzyme_rate_equation = strip_time(CPS.enzyme_rate_error(False)) print("Enzyme Net Flux Equation\n" + "-"*24) sym.pprint(enzyme_rate_equation) # Solve for last unknown concentration symbolically sol = sym.solveset(enzyme_rate_equation.subs(enzyme_solutions), "cps_m") # Update solution dictionary with the new solution enzyme_solutions[sym.Symbol("cps_m")] = list(sol)[0] # Update solutions with free variable solutions enzyme_solutions = { enzyme_module_form: sym.simplify(solution.subs(enzyme_solutions)) for enzyme_module_form, solution in enzyme_solutions.items()} args = set() for sol in enzyme_solutions.values(): args.update(sol.atoms(sym.Symbol)) print("\n", args)

numerical_values = {} # Get ligand IDs and parameter IDs ligand_ids = sorted([str(ligand) for ligand in CPS.enzyme_module_ligands]) parameter_ids = ["_".join((CPS.id, abbrev)) for abbrev in abbreviations + ["L"]] print("Ligand IDs: " + str(ligand_ids)) print("Parameter IDs: " + str(parameter_ids)) # Get the steady state flux value and add to numerical values CPS.enzyme_rate = model.reactions.get_by_id(CPS.id).steady_state_flux numerical_values.update({CPS.enzyme_flux_symbol_str: CPS.enzyme_rate}) # Get the steady state concentration values and add to numerical values initial_conditions = { str(ligand): model.initial_conditions[model.metabolites.get_by_id(ligand)] for ligand in ligand_ids} # Define parameter values and add to numerical values # Because of the QEA, invert dissociation constants for Keq parameter_values = { "Keq_" + parameter_id: value for parameter_id, value in zip(parameter_ids, [1/1, 1/4, 1/0.5, 0.0011])} # Display numerical values print("\nNumerical Values\n----------------") for k, v in numerical_values.items(): print("{0} = {1}".format(k, v))

# Match abbreviations to their corresponding ligands abbreviation_dict = {"CPS_N": "nh4_m", "CPS_H": "hco3_m", "CPS_A": "atp", "CPS_L": ""} k2K = {sym.Symbol("kr_" + p): sym.Symbol("kf_" + p)*sym.Symbol("K_" + p) for p in abbreviation_dict.keys()} enzyme_solutions = {met: sym.simplify(Keq2k(solution).subs(enzyme_solutions).subs(k2K)) for met, solution in enzyme_solutions.items()} K_values = dict(zip(["K_" + p for p in abbreviation_dict], [1, 4, 0.5, 0.0011])) for abbrev, ligand_id in abbreviation_dict.items(): K_str = "K_" + abbrev if ligand_id: numerical_value = K_values[K_str]/initial_conditions[ligand_id] else: numerical_value = 1/K_values[K_str] numerical_values[sym.Symbol(K_str)] = numerical_value enzyme_solutions = {met: sym.simplify(solution.subs(numerical_values)) for met, solution in enzyme_solutions.items()} # Display numerical values print("\nNumerical Values\n----------------") for k, v in numerical_values.items(): print("{0} = {1}".format(k, v))

# Substitute values into equations enzyme_solutions = { enzyme_module_form: sym.simplify(solution.subs(numerical_values)) for enzyme_module_form, solution in enzyme_solutions.items()} args = set() for sol in enzyme_solutions.values(): args.update(sol.atoms(sym.Symbol)) print(args)

sym.pprint(CPS.enzyme_concentration_total_equation)

CPS.enzyme_concentration_total = 40 print(CPS.enzyme_concentration_total)

enzyme_total_constraint = abs(strip_time(CPS.enzyme_concentration_total_error(use_values=False))) sym.pprint(enzyme_total_constraint)

# Substitute value for enzyme concentration total enzyme_total_constraint = enzyme_total_constraint.subs({CPS.enzyme_total_symbol_str: CPS.enzyme_concentration_total}) # Substitute solutions into constraint and simplify enzyme_total_constraint = sym.simplify(enzyme_total_constraint.subs(enzyme_solutions)) sym.pprint(enzyme_total_constraint)

# Create a sorted tuple of the arguments to ensure the input format does not change args = tuple(sorted([str(arg) for arg in list(args)])) # Create the objective function as a lambda function objective_function = lambda x: sym.lambdify(args, enzyme_total_constraint)(*x)

# Set the values for the constraint bounds r_T_lb, r_T_ub = (0.01, 0.015) # Make a symbolic expression for enzyme fraction. r_T_expr = CPS.make_enzyme_fraction( categorized_attr="forms", top="Relaxed", bottom="Equation", use_values=False) # Substitute solutions into the expression to make # solely dependent on the rate constants r_T_expr = sym.simplify(strip_time(r_T_expr).subs(enzyme_solutions)) # Make lambda functions for the T fraction constraint r_T_lb_constraint = lambda x: sym.lambdify(args, r_T_expr - r_T_lb)(*x) r_T_ub_constraint = lambda x: sym.lambdify(args, r_T_ub - r_T_expr)(*x)

print("Ordered Args: {0}\n".format(str(args))) # Set arguments for minimization kf_bounds = ((1e2, 1e8), (1e2, 1e8), (1e2, 1e8), (1e2, 1e8)) initial_guess = [ 3.07e5, 2e5, 1e6, 2e6] # Find a feasible solution sol = optimize.minimize( objective_function, x0=initial_guess, method="trust-constr", bounds=kf_bounds, options={"gtol": 1e-20, "xtol": 1e-20, "maxiter": 1e4, "disp": True}) # Check whether optimzation was successful print("\nOptimization Success: {0}".format(sol.success)) if sol.success: # Update the paramter values dictionary with the feasible solution parameter_values.update(dict(zip(args, [round(x) for x in sol.x]))) print("Optimization Optimality: {0:.4e}".format(sol.optimality)) print("Parameter Solutions: {:}".format(str({arg: parameter_values[arg] for arg in args}))) # Plug solutions back into constraints for validation print("Optimization Error: {0:.4e}".format(enzyme_total_constraint.subs(parameter_values)))

# Display numerical values for k, v in parameter_values.items(): print("{0} = {1}".format(k, v))

# Substitute values into equations initial_conditions.update({ str(enzyme_module_form): float(sym.simplify(solution.subs(parameter_values))) for enzyme_module_form, solution in enzyme_solutions.items()}) for header, dictlist in zip(["Ligand", "\nEnzyme"], [CPS.enzyme_module_ligands, CPS.enzyme_module_forms]): header += " Concentrations" print("\n".join([header, "-" * len(header)])) for form in dictlist: ic = initial_conditions[form.id] print("{0} = {1}".format(form.id, ic))

# Set initial conditions for met, concentration in initial_conditions.items(): CPS.metabolites.get_by_id(str(met)).ic = concentration # Add the custom parameters and values for kf and Keq to model CPS.custom_parameters.update(parameter_values) # CPS_L = CPS.enzyme_module_reactions.CPS_L # CPS_L.kf = CPS.custom_parameters.pop(CPS_L.kf_str) # CPS_L.Keq = CPS.custom_parameters.pop(CPS_L.Keq_str) # Set parameter values in reaction fields for group in CPS.enzyme_module_reactions_categorized: if group.id == "nh4_m_binding": param_id = "CPS_N" elif group.id == "hco3_m_binding": param_id = "CPS_H" elif group.id == "atp_binding": param_id = "CPS_A" elif group.id == "catalyzation": param_id = "CPS" else: continue for reaction in group.members: kf, Keq = ("kf_" + param_id, "Keq_" + param_id) if kf in CPS.custom_parameters: reaction.kf = CPS.custom_parameters[kf] if Keq in CPS.custom_parameters: reaction.Keq = CPS.custom_parameters[Keq]

CPS.update_S(array_type="DataFrame", dtype=int)

qcqa_model(CPS, parameters=True, concentrations=True, fluxes=False, superfluous=True, elemental=True)

t_fraction = CPS.make_enzyme_fraction("forms", top="Relaxed", bottom="Equation", use_values=True) print("Enzyme T-fraction: {:.4f}".format(t_fraction)) print("Concentration Absolute Error: {0:.4e}".format( abs(CPS.enzyme_concentration_total_error(use_values=True)))) print("Flux Absolute Error: {0:.4e}".format( abs(CPS.enzyme_rate_error(use_values=True))))

# Load and merge main model with CPS module # Remove the CPS MassReaction, then merge the EnzymeModule into the MassModel model.remove_reactions([model.reactions.get_by_id("CPS")]) model_CPS = model.merge(CPS) model_CPS

print(str(model_CPS.enzyme_modules) + "\n") print("Attribute Accessors:\n-------------------\n" + "\n".join(list( model_CPS.enzyme_modules.CPS.keys())) + "\n") model_CPS.enzyme_modules.CPS

# Setup simulation object, ensure model is at steady state sim = Simulation(model_CPS, verbose=True) sim.find_steady_state(model_CPS, strategy="simulate", update_values=True, verbose=True, tfinal=1e4, steps=1e6) # Simulate from 0 to 1000 with 10001 points in the output conc_sol, flux_sol = sim.simulate( model_CPS,time=(0, 1e3, 1e4 + 1)) # Quickly render and display time profiles conc_sol.view_time_profile()

# Plot time profile of ammonium and urea plot_time_profile( conc_sol, observable=["nh4_m", "carb_phos_m", "ornithine", "citrulline", "aspartate_c", "argi_suc_c", "urea_c"], legend="right outside", plot_function="semilogx", xlabel="Time", ylabel="Concentration", title=("Concentration Solutions", {"size": "large"}))

# Plot time profil for fluxes on log scale plot_time_profile( flux_sol, legend="right outside", plot_function="semilogx", xlabel="Time [hr]", ylabel="Flux [mM/hr]", title=("Flux Profile", L_FONT));