<jupyter_start><jupyter_text>Getting Started with MASSpyIn this notebook example, objects essential to **MASSpy** are explored.  ModelsIn **MASSpy**, a model is represented by a ``mass.MassModel`` object. **MASSpy** comes bundled with example models, including a "textbook" model$^1$ of human red blood cell metabolism. To load this example model:<jupyter_code>from operator import attrgetter

import mass
import mass.example_data

model = mass.example_data.create_example_model("textbook")<jupyter_output>Set parameter Username<jupyter_text>Several attributes of the `MassModel`, including model reactions and metabolites, are special types of lists that contain objects related to the model. Each specialized list is called a `cobra.DictList` and is made up of the corresponding objects. For example, the `reactions` attribute contains the `mass.MassReaction` objects, and the `metabolites` attribute contains the `mass.MassMetabolite` objects.<jupyter_code>print("Number of metabolites: " + str(len(model.metabolites)))
print("Number of reactions: " + str(len(model.reactions)))<jupyter_output>Number of metabolites: 68
Number of reactions: 76<jupyter_text>When using a [Jupyter notebook](https://jupyter-notebook-beginner-guide.readthedocs.io/en/latest/), this type of information is rendered as a table.<jupyter_code>model<jupyter_output><empty_output><jupyter_text>Just like a regular list, objects in a `cobra.DictList` can be retrieved by index. For example, to get the 30th reaction in the model (at index 29 because of [0-indexing](https://en.wikipedia.org/wiki/Zero-based_numbering)):<jupyter_code>model.reactions[29]<jupyter_output><empty_output><jupyter_text>Items also can be retrieved by their `id` attribute using the `cobra.DictList.get_by_id()` method. For example, to get the cytosolic atp metabolite object with identifier "atp_c":<jupyter_code>model.metabolites.get_by_id("atp_c")<jupyter_output><empty_output><jupyter_text>If care is taken when assigning object identifiers (e.g., does not start with a number, does not contain certain characters such as "-"), it is possible to access objects inside of a `cobra.DictList` as if they were attributes.<jupyter_code>print(model.reactions.DPGase)<jupyter_output>DPGase: _23dpg_c + h2o_c --> _3pg_c + pi_c<jupyter_text>To ensure all identifiers comply with [Systems Biology Markup Language (SBML)](http://sbml.org/Main_Page) and allow for interactive use, utilizing the identifiers from the [BiGG Models database](http://bigg.ucsd.edu/) is highly recommended.Guidelines for BiGG identifiers are found [here](https://github.com/SBRG/bigg_models/wiki/BiGG-Models-ID-Specification-and-Guidelines).  ReactionsIn **MASSpy**, a reaction is represented by a ``mass.MassReaction`` object. A particular reaction can be retrieved by its `id` using the `cobra.DictList.get_by_id()` method. Below, the reaction with identifier "PGI" is inspected.<jupyter_code>PGI = model.reactions.get_by_id("PGI")
PGI<jupyter_output><empty_output><jupyter_text>The full name and the chemical reaction are viewed as strings. If defined, the flux value for the reaction at steady state also can be viewed.<jupyter_code>print(PGI.name)
print(PGI.reaction)
print(PGI.steady_state_flux)<jupyter_output>Glucose-6-phosphate isomerase
g6p_c <=> f6p_c
0.9098871145647632<jupyter_text>The symbolic rate equation for the reaction is viewed using the `rate` attribute. The rate is returned as a **SymPy** symbolic expression.<jupyter_code>print(PGI.rate)<jupyter_output>kf_PGI*(g6p_c(t) - f6p_c(t)/Keq_PGI)<jupyter_text>The above rate is considered a Type 1 rate law because of the reaction parameters used. There are three types of mass action rate laws that can be generated;  * Type 1 rates utilize the forward rate and equilibrium constants.  * Type 2 rates utilize the forward rate and reverse rate constants.  * Type 3 rates utilize the equilibrium and reverse rate constants.  To view the reaction rate as a Type 2 rate equation, the `get_mass_action_rate()` method can be used.<jupyter_code>print(PGI.get_mass_action_rate(rate_type=2))<jupyter_output>kf_PGI*g6p_c(t) - kr_PGI*f6p_c(t)<jupyter_text>Because the reaction has its `reversible` attribute set as `True`, the net mass action rate contains both forward and reverse rate components.<jupyter_code>print("Reversible: {0!r}".format(PGI.reversible))
print("Forward rate: {0!r}".format(
    PGI.get_forward_mass_action_rate_expression(rate_type=2)))
print("Reverse rate: {0!r}".format(
    PGI.get_reverse_mass_action_rate_expression(rate_type=2)))
print("Net reaction rate: {0!r}".format(
    PGI.get_mass_action_rate(rate_type=2)))<jupyter_output>Reversible: True
Forward rate: kf_PGI*g6p_c(t)
Reverse rate: kr_PGI*f6p_c(t)
Net reaction rate: kf_PGI*g6p_c(t) - kr_PGI*f6p_c(t)<jupyter_text>To view the defined parameters for the rate and equilibrium constants, the `parameters` attribute is used to return a `dict` containing the parameter identifiers and their values.<jupyter_code>PGI.parameters<jupyter_output><empty_output><jupyter_text>Parameter identifiers for a reaction can be obtained using various attributes.<jupyter_code>print(PGI.flux_symbol_str)
print(PGI.kf_str)
print(PGI.Keq_str)
print(PGI.kr_str)<jupyter_output>v_PGI
kf_PGI
Keq_PGI
kr_PGI<jupyter_text>Changing the `reversible` attribute of the reaction affects the net rate equation:<jupyter_code>PGI.reversible = False
print("Reversible: {0!r}".format(PGI.reversible))
print("Net reaction rate: {0!r}".format(PGI.rate))<jupyter_output>Reversible: False
Net reaction rate: kf_PGI*g6p_c(t)<jupyter_text>Note that changing the `reversible` attribute of the reaction can affect the parameters:<jupyter_code>PGI.parameters<jupyter_output><empty_output><jupyter_text>The reaction can be checked for whether it is mass balanced using the `check_mass_balance()` method. This method returns the elements that violate mass balance. If it comes back empty, then the reaction is mass balanced.<jupyter_code>PGI.check_mass_balance()<jupyter_output><empty_output><jupyter_text>The `add_metabolites()` method can be used to add metabolites to a reaction by passing in a `dict` that contains the `MassMetabolite` objects and their coefficients:<jupyter_code>PGI.add_metabolites({model.metabolites.get_by_id("h_c"): -1})
print(PGI)<jupyter_output>PGI: g6p_c + h_c --> f6p_c<jupyter_text>The reaction is no longer mass balanced.<jupyter_code>PGI.subtract_metabolites({model.metabolites.get_by_id("h_c"): -1})
print(PGI)
print(PGI.check_mass_balance())<jupyter_output>PGI: g6p_c --> f6p_c
{}<jupyter_text>MetabolitesIn **MASSpy**, a metabolite is represented by a ``MassMetabolite`` object. A particular metabolite can be retrieved by its `id` using the `cobra.DictList.get_by_id()` method. Below, the cytosolic glucose 6-phosphate metabolite with identifier "g6p_c" is inspected.<jupyter_code>g6p_c = model.metabolites.get_by_id("g6p_c")
g6p_c<jupyter_output><empty_output><jupyter_text>The full name and the compartment where the metabolite is located ("c" for cytosol) are viewed as strings:<jupyter_code>print(g6p_c.name)
print(g6p_c.compartment)<jupyter_output>D-Glucose 6-phosphate
c<jupyter_text>The chemical formula and associated charge of the metabolite also can be viewed:<jupyter_code>print(g6p_c.formula)
print(g6p_c.charge)<jupyter_output>C6H11O9P
-2<jupyter_text>Reactions in which the metabolite participates are obtained as a `frozenset` from the `reactions` attribute. This can be used to count the number of reactions that utilize the metabolite.<jupyter_code>print("Number of reactions involving {0}: {1}".format(
    g6p_c.id, len(g6p_c.reactions)))<jupyter_output>Number of reactions involving g6p_c: 3<jupyter_text>The ordinary differential equation (ODE), which represents the change in metabolite concentration over time, is determined by reactions that consume or produce the metabolite. The `oridinary_differential_equation` attribute is used to view the current ODE for the metabolite. The ODE is returned as a symbolic expression.<jupyter_code>print(g6p_c.ordinary_differential_equation)<jupyter_output>-kf_G6PDH2r*(g6p_c(t)*nadp_c(t) - _6pgl_c(t)*nadph_c(t)/Keq_G6PDH2r) + kf_HEX1*(atp_c(t)*glc__D_c(t) - adp_c(t)*g6p_c(t)/Keq_HEX1) - kf_PGI*g6p_c(t)<jupyter_text>Numerical solutions are obtained by integrating ODEs. To integrate an ODE, a metabolite concentration at time $t = 0$ must be defined as an initial condition. Initial conditions are accessed and changed using the `initial_condition` attribute.<jupyter_code>g6p_c.initial_condition = 0.8
print(g6p_c.initial_condition)<jupyter_output>0.8<jupyter_text>Certain attributes have alias attribute accessors. For example, the `initial_condition` and `ordinary_differential_equation` attributes can be accessed via `ic` and `ode`, respectively.<jupyter_code>print(g6p_c.ic)<jupyter_output>0.8<jupyter_text>The `fixed` attribute indicates whether the concentration of the metabolite is allowed to vary over time, with `True` meaning that the metabolite's initial condition is treated as a constant concentration value. Fixed metabolites have an ODE equal to 0.<jupyter_code>g6p_c.fixed = True
print(g6p_c.ode)<jupyter_output>0<jupyter_text>Additional Model ObjectsThe following are additional objects that are stored within a `MassModel`. Unlike metabolites and reactions, which are essential for defining the system of ODEs, these objects are not always necessary for dynamic simulation of models.However, these objects are still important to **MASSpy** and have a variety of uses, which include: aiding in the management of large models, sharing models among users, tracking additional information relevant to the system, and enabling various workflows for genome-scale kinetic models.  GenesBecause the `mass.MassReaction` inherits from the `cobra.Reaction`, **MASSpy** is also capable of handling genes and gene-protein-reaction (GPR) relationships. Note that **MASSpy** directly utilizes the `cobra.Gene` object for the representation and management of genes.The `gene_reaction_rule` is a Boolean representation of the gene requirements for this reaction to be active (Schellenberger et. al, 2011).GPRs are stored as the `gene_reaction_rule` of reaction objects.Altering a `gene_reaction_rule` will create new gene objects, if necessary, and update all relationships.<jupyter_code>PGI.gene_reaction_rule = "New_Gene"
PGI.gene_reaction_rule<jupyter_output><empty_output><jupyter_text>`Gene` objects are returned from a reaction as a `frozenset`:<jupyter_code>PGI.genes<jupyter_output><empty_output><jupyter_text>Newly created genes are added to the model upon creation. The gene objects are stored in a `cobra.DictList` as a part of the `genes` attribute. To access a gene from the model:<jupyter_code>new_gene = model.genes.get_by_id("New_Gene")
new_gene<jupyter_output><empty_output><jupyter_text>Gene objects are tracked by both its associated reaction objects and the the model. Changing a reaction's `gene_reaction_rule` may remove the gene's association from the reaction, but it does not remove the gene from the model.<jupyter_code>PGI.gene_reaction_rule = ""
print("Reaction Genes: {0!r}".format(PGI.genes))
new_gene<jupyter_output>Reaction Genes: frozenset()<jupyter_text>EnzymeModulesA `mass.EnzymeModule` is a specialized `MassModel` that represents a reconstruction of an enzyme's mechanism. Upon merging an `EnzymeModule` into a `MassModel`, the `EnzymeModule` is converted into a `mass.EnzymeModuleDict`, a specialized dictionary object. The `EnzymeModuleDict` is subsequently stored in a `cobra.DictList` and is accessible through the `enzyme_modules` attribute.<jupyter_code>PFK = model.enzyme_modules.get_by_id("PFK")
PFK<jupyter_output><empty_output><jupyter_text>The process of creating an `EnzymeModuleDict` preserves information stored in various attributes specific to the enzyme module and allows them to be accessed quickly after the merging process. Because the `EnzymeModuleDict` inherits from an `OrderedDict`, it has the same methods and behaviors.<jupyter_code>for key, value in PFK.items():
    if isinstance(value, list):
        print("{0}: {1}".format(key, len(value)))<jupyter_output>enzyme_module_ligands: 6
enzyme_module_forms: 20
enzyme_module_reactions: 24
enzyme_module_ligands_categorized: 5
enzyme_module_forms_categorized: 5
enzyme_module_reactions_categorized: 6<jupyter_text>`EnzymeModuleDict` objects also can have their contents accessed by using its `dict` keys as attribute accessors:<jupyter_code>print(PFK.id)
print(PFK.subsystem)
print(PFK.enzyme_concentration_total)<jupyter_output>PFK
Glycolysis
3.3e-05<jupyter_text>See the section on [EnzymeModules](./enzyme_modules.ipynb) for more information on working with `EnzymeModule` and related objects.  GroupsGroups are objects for holding information regarding pathways, subsystems, or any custom grouping of objects within a `MassModel`. **MASSpy** directly utilizes the `cobra.Group` object, which are implemented based on the [SBML Group](http://sbml.org/Software/libSBML/5.18.0/docs/python-api/classlibsbml_1_1_group.html) specifications.`Group` objects are stored in a `cobra.DictList` as the `groups` attribute of the `MassModel`.<jupyter_code>print("Number of groups: {0}".format(len(model.groups)))<jupyter_output>Number of groups: 16<jupyter_text>There are several different ways to work with group objects. One way that **MASSpy** utilizes group objects is to aid with categorizing and grouping various objects associated with `EnzymeModule` objects.For example, the group `Products` contains all of metabolites that are the products of the reaction catalyzed by the `PFK` enzyme. A `set` containing the associated metabolite objects is returned by the `members` attribute.<jupyter_code>products = model.groups.get_by_id("Products")
products.members<jupyter_output><empty_output><jupyter_text>Groups are also used to categorize reactions. For example, the group `atp_c_binding` contains all of the reactions that represent the binding of ATP to the free active sites of the PFK enzyme.<jupyter_code>complexed_w_atp = model.groups.get_by_id("atp_c_binding")
for member in complexed_w_atp.members:
    print(member)<jupyter_output>PFK_R41: atp_c + pfk_R4_c <=> pfk_R4_A_c
PFK_R01: atp_c + pfk_R0_c <=> pfk_R0_A_c
PFK_R11: atp_c + pfk_R1_c <=> pfk_R1_A_c
PFK_R31: atp_c + pfk_R3_c <=> pfk_R3_A_c
PFK_R21: atp_c + pfk_R2_c <=> pfk_R2_A_c<jupyter_text>Because groups are sets, members of groups are returned in no particular order. To maintain a consistent order, the `sorted()` function is used with the `attrgetter()` function from the `operator` module to sort members by a particular object attribute.<jupyter_code>for member in sorted(complexed_w_atp.members, key=attrgetter("id")):
    print(member)<jupyter_output>PFK_R01: atp_c + pfk_R0_c <=> pfk_R0_A_c
PFK_R11: atp_c + pfk_R1_c <=> pfk_R1_A_c
PFK_R21: atp_c + pfk_R2_c <=> pfk_R2_A_c
PFK_R31: atp_c + pfk_R3_c <=> pfk_R3_A_c
PFK_R41: atp_c + pfk_R4_c <=> pfk_R4_A_c<jupyter_text>Units`Unit` and `UnitDefinition` objects are implemented, as per the [SBML Unit](http://sbml.org/Software/libSBML/5.18.0/docs/python-api/classlibsbml_1_1_unit.html) and [SBML UnitDefinition](http://sbml.org/Software/libSBML/5.18.0/docs/python-api/classlibsbml_1_1_unit_definition.html) specifications. The primary purpose of these objects is to inform users of the model's units, providing context to model values and observed results.It is important to note that unit consistency is **NOT** checked by the `MassModel`, meaning that it is incumbent upon users to maintain consistency of units and associated numerical values in a model. `UnitDefinition` objects are stored in a `cobra.DictList`, accessible through the `units` attribute.<jupyter_code>model.units<jupyter_output><empty_output><jupyter_text>`UnitDefinition` objects have identifiers and optional names. Therefore, a specific unit can be accessed using the `get_by_id()` method.<jupyter_code>concentration_unit = model.units.get_by_id("mM")
print(concentration_unit.id)
print(concentration_unit.name)<jupyter_output>mM
Millimolar<jupyter_text>A `UnitDefinition` is comprised of base units, which are stored in the `list_of_units` attribute. Each unit must have a defined `kind`, `exponent`, `scale`, and `multiplier`.<jupyter_code>for unit in concentration_unit.list_of_units:
    print(unit)<jupyter_output>kind: litre; exponent: -1; scale: 0; multiplier: 1
kind: mole; exponent: 1; scale: -3; multiplier: 1<jupyter_text>**MASSpy** contains some commonly defined `Unit` objects to aid in the creation of `UnitDefinition` objects, which are viewed using the `print_defined_unit_values()` function.<jupyter_code>mass.core.units.print_defined_unit_values("Units")<jupyter_output>╒════════════════════════════════════════════════════════════════════╕
│ Pre-defined Units                                                  │
╞════════════════════════════════════════════════════════════════════╡
│ Unit        Definition                                             │
│ ----------  ------------------------------------------------------ │
│ mole        kind: mole; exponent: 1; scale: 0; multiplier: 1       │
│ millimole   kind: mole; exponent: 1; scale: -3; multiplier: 1      │
│ litre       kind: litre; exponent: 1; scale: 0; multiplier: 1      │
│ per_litre   kind: litre; exponent: -1; scale: 0; multiplier: 1     │
│ second      kind: second; exponent: 1; scale: 0; multiplier: 1     │
│ per_second  kind: second; exponent: -1; scale: 0; multiplier: 1    │
│ hour        kind: second; exponent: 1; scale: 0; multiplier: 3600  │
│ per_hour    kind: second; exponent: -1; scale: 0; multiplier: 3600 │
│ per_gDW     kind: gram; exponent: -1; scale: 0; multiplier: 1      │
╘═════[...]
{'fraction_non_alphanumeric': 0.10099300383660573, 'fraction_numerical': 0.013033175355450236, 'mean_word_length': 6.2792607802874745, 'pattern_counts': {'":': 1, '<': 149, '<?xml version=': 0, '>': 152, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '41', 'dir': 'jupyter-structured-clean-dedup', 'id': '256420', 'license': 'MIT', 'n_tokens_mistral': 6036, 'n_tokens_neox': 5642, 'n_words': 1763, 'path': 'docs/tutorials/getting_started_with_masspy.ipynb', 'repo_name': 'z-haiman/MASSpy', 'size': '17724'}