"""Define the Model class."""
import logging
from copy import copy, deepcopy
from functools import partial
from types import ModuleType
from typing import TYPE_CHECKING, Dict, Iterable, List, Optional, Tuple, Union
from warnings import warn
import optlang
from optlang.symbolics import Basic, Zero
from ..medium import find_boundary_types, find_external_compartment, sbo_terms
from ..util.context import HistoryManager, get_context, resettable
from ..util.solver import (
add_cons_vars_to_problem,
assert_optimal,
check_solver,
interface_to_str,
remove_cons_vars_from_problem,
set_objective,
)
from ..util.util import AutoVivification, format_long_string
from .configuration import Configuration
from .dictlist import DictList
from .gene import Gene
from .group import Group
from .metabolite import Metabolite
from .object import Object
from .reaction import Reaction
from .solution import get_solution
if TYPE_CHECKING:
import pandas as pd
from optlang.container import Container
from cobra import Solution
from cobra.summary import ModelSummary
from cobra.util.solver import CONS_VARS
[docs]logger = logging.getLogger(__name__)
[docs]configuration = Configuration()
[docs]class Model(Object):
"""Class representation for a cobra model.
Parameters
----------
id_or_model: str or Model, optional
String to use as model id, or actual model to base new model one.
If string, it is used as id. If model, a new model object is
instantiated with the same properties as the original model (default None).
name: str, optional
Human readable string to be model description (default None).
Attributes
----------
reactions : DictList
A DictList where the key is the reaction identifier and the value a
Reaction
metabolites : DictList
A DictList where the key is the metabolite identifier and the value a
Metabolite
genes : DictList
A DictList where the key is the gene identifier and the value a
Gene
groups : DictList
A DictList where the key is the group identifier and the value a
Group
"""
def __init__(
self, id_or_model: Union[str, "Model", None] = None, name: Optional[str] = None
) -> None:
"""Initialize the Model."""
if isinstance(id_or_model, Model):
Object.__init__(self, name=name)
self.__setstate__(id_or_model.__dict__)
self._solver = id_or_model.solver
else:
Object.__init__(self, id_or_model, name=name)
self.genes = DictList()
self.reactions = DictList() # A list of cobra.Reactions
self.metabolites = DictList() # A list of cobra.Metabolites
self.groups = DictList() # A list of cobra.Groups
# genes based on their ids {Gene.id: Gene}
self._compartments = {}
self._contexts = []
# from cameo ...
# if not hasattr(self, '_solver'): # backwards compatibility
# with older cobrapy pickles?
interface = check_solver(configuration.solver)
self._solver = interface.Model()
self._solver.objective = interface.Objective(Zero)
self._populate_solver(self.reactions, self.metabolites)
self._tolerance = None
self.tolerance = configuration.tolerance
[docs] def __setstate__(self, state: Dict) -> None:
"""Make sure all cobra.Objects in the model point to the model.
Parameters
----------
state: dict
"""
self.__dict__.update(state)
for y in ["reactions", "genes", "metabolites"]:
for x in getattr(self, y):
x._model = self
if not hasattr(self, "name"):
self.name = None
[docs] def __getstate__(self) -> Dict:
"""Get state for serialization.
Ensures that the context stack is cleared prior to serialization,
since partial functions cannot be pickled reliably.
Returns
-------
odict: Dict
A dictionary of state, based on self.__dict__.
"""
odict = self.__dict__.copy()
odict["_contexts"] = []
return odict
@property
[docs] def solver(self) -> "optlang.interface.Model":
"""Get the attached solver instance.
The associated the solver object, which manages the interaction with
the associated solver, e.g. glpk.
This property is useful for accessing the optimization problem
directly and to define additional non-metabolic constraints.
Examples
--------
>>> from cobra.io import load_model
>>> model = load_model("textbook")
>>> new = model.problem.Constraint(model.objective.expression, lb=0.99)
>>> model.solver.add(new)
"""
return self._solver
@solver.setter
@resettable
def solver(self, value: Union[str, ModuleType]) -> None:
"""Set the attached solver instance.
The associated the solver object, which manages the interaction with
the associated solver, e.g. glpk.
This property is useful for accessing the optimization problem
directly and to define additional non-metabolic constraints.
Parameters
----------
value: ModuleType or str
The solver to set, which will be checked by `check_solver()`.
"""
interface = check_solver(value)
# Do nothing if the solver did not change
if self.problem == interface:
return
self._solver = interface.Model.clone(self._solver)
@property
[docs] def tolerance(self) -> float:
"""Get the tolerance.
Returns
-------
float
The tolerance of the mdoel.
"""
return self._tolerance
@tolerance.setter
def tolerance(self, value: float) -> None:
"""Set the tolerance.
Parameters
----------
value: float
Value to set tolerance.
"""
solver_tolerances = self.solver.configuration.tolerances
try:
solver_tolerances.feasibility = value
except AttributeError:
logger.info(
f"The current solver interface {interface_to_str(self.problem)} "
f"doesn't support setting the feasibility tolerance."
)
try:
solver_tolerances.optimality = value
except AttributeError:
logger.info(
f"The current solver interface {interface_to_str(self.problem)} "
f"doesn't support setting the optimality tolerance."
)
try:
solver_tolerances.integrality = value
except AttributeError:
logger.info(
f"The current solver interface {interface_to_str(self.problem)} "
f"doesn't support setting the integrality tolerance."
)
self._tolerance = value
@property
[docs] def compartments(self) -> Dict:
"""Return all metabolites' compartments.
Returns
-------
dict
A dictionary of metabolite compartments, where the keys are the short
version (one letter version) of the compartmetns, and the values are the
full names (if they exist).
"""
return {
met.compartment: self._compartments.get(met.compartment, "")
for met in self.metabolites
if met.compartment is not None
}
@compartments.setter
def compartments(self, value: Dict) -> None:
"""Get or set the dictionary of current compartment descriptions.
Assigning a dictionary to this property updates the model's
dictionary of compartment descriptions with the new values.
Parameters
----------
value : dict
Dictionary mapping compartments abbreviations to full names.
Examples
--------
>>> from cobra.io import load_model
>>> model = load_model("textbook")
>>> model.compartments = {'c': 'the cytosol'}
>>> model.compartments
{'c': 'the cytosol', 'e': 'extracellular'}
"""
self._compartments.update(value)
@property
[docs] def medium(self) -> Dict[str, float]:
"""Get the constraints on the model exchanges.
`model.medium` returns a dictionary of the bounds for each of the
boundary reactions, in the form of `{rxn_id: bound}`, where `bound`
specifies the absolute value of the bound in direction of metabolite
creation (i.e., lower_bound for `met <--`, upper_bound for `met -->`)
Returns
-------
Dict[str, float]
A dictionary with rxn.id (str) as key, bound (float) as value.
"""
def is_active(reaction: Reaction) -> bool:
"""Determine if boundary reaction permits flux towards creating metabolites.
Parameters
----------
reaction: cobra.Reaction
Returns
-------
bool
True if reaction produces metaoblites and has upper_bound above 0
or if reaction consumes metabolites and has lower_bound below 0 (so
could be reversed).
"""
return (bool(reaction.products) and (reaction.upper_bound > 0)) or (
bool(reaction.reactants) and (reaction.lower_bound < 0)
)
def get_active_bound(reaction: Reaction) -> float:
"""For an active boundary reaction, return the relevant bound.
Parameters
----------
reaction: cobra.Reaction
Returns
-------
float:
upper or minus lower bound, depenending if the reaction produces or
consumes metaoblties.
"""
if reaction.reactants:
return -reaction.lower_bound
elif reaction.products:
return reaction.upper_bound
return {
rxn.id: get_active_bound(rxn) for rxn in self.exchanges if is_active(rxn)
}
@medium.setter
def medium(self, medium: Dict[str, float]) -> None:
"""Set the constraints on the model exchanges.
`model.medium` returns a dictionary of the bounds for each of the
boundary reactions, in the form of `{rxn_id: rxn_bound}`, where `rxn_bound`
specifies the absolute value of the bound in direction of metabolite
creation (i.e., lower_bound for `met <--`, upper_bound for `met -->`)
Parameters
----------
medium: dict
The medium to initialize. medium should be a dictionary defining
`{rxn_id: bound}` pairs.
"""
def set_active_bound(reaction: Reaction, bound: float) -> None:
"""Set active bound.
Parameters
----------
reaction: cobra.Reaction
Reaction to set
bound: float
Value to set bound to. The bound is reversed and set as lower bound
if reaction has reactants (metabolites that are consumed). If reaction
has reactants, it seems the upper bound won't be set.
"""
if reaction.reactants:
reaction.lower_bound = -bound
elif reaction.products:
reaction.upper_bound = bound
# Set the given media bounds
media_rxns = []
exchange_rxns = frozenset(self.exchanges)
for rxn_id, rxn_bound in medium.items():
rxn = self.reactions.get_by_id(rxn_id)
if rxn not in exchange_rxns:
logger.warning(
f"{rxn.id} does not seem to be an an exchange reaction. "
f"Applying bounds anyway."
)
media_rxns.append(rxn)
# noinspection PyTypeChecker
set_active_bound(rxn, rxn_bound)
frozen_media_rxns = frozenset(media_rxns)
# Turn off reactions not present in media
for rxn in exchange_rxns - frozen_media_rxns:
is_export = rxn.reactants and not rxn.products
set_active_bound(
rxn, min(0.0, -rxn.lower_bound if is_export else rxn.upper_bound)
)
[docs] def copy(self) -> "Model":
"""Provide a partial 'deepcopy' of the Model.
All the Metabolite, Gene, and Reaction objects are created anew but
in a faster fashion than deepcopy.
Returns
-------
cobra.Model: new model copy
"""
new = self.__class__()
do_not_copy_by_ref = {
"metabolites",
"reactions",
"genes",
"notes",
"annotation",
"groups",
}
for attr in self.__dict__:
if attr not in do_not_copy_by_ref:
new.__dict__[attr] = self.__dict__[attr]
new.notes = deepcopy(self.notes)
new.annotation = deepcopy(self.annotation)
new.metabolites = DictList()
do_not_copy_by_ref = {"_reaction", "_model"}
for metabolite in self.metabolites:
new_met = metabolite.__class__()
for attr, value in metabolite.__dict__.items():
if attr not in do_not_copy_by_ref:
new_met.__dict__[attr] = copy(value) if attr == "formula" else value
new_met._model = new
new.metabolites.append(new_met)
new.genes = DictList()
for gene in self.genes:
new_gene = gene.__class__(None)
for attr, value in gene.__dict__.items():
if attr not in do_not_copy_by_ref:
new_gene.__dict__[attr] = (
copy(value) if attr == "formula" else value
)
new_gene._model = new
new.genes.append(new_gene)
new.reactions = DictList()
do_not_copy_by_ref = {"_model", "_metabolites", "_genes"}
for reaction in self.reactions:
new_reaction = reaction.__class__()
for attr, value in reaction.__dict__.items():
if attr not in do_not_copy_by_ref:
new_reaction.__dict__[attr] = copy(value)
new_reaction._model = new
new.reactions.append(new_reaction)
# update awareness
for metabolite, stoic in reaction._metabolites.items():
new_met = new.metabolites.get_by_id(metabolite.id)
new_reaction._metabolites[new_met] = stoic
new_met._reaction.add(new_reaction)
new_reaction.update_genes_from_gpr()
new.groups = DictList()
do_not_copy_by_ref = {"_model", "_members"}
# Groups can be members of other groups. We initialize them first and
# then update their members.
for group in self.groups:
new_group: Group = group.__class__(group.id)
for attr, value in group.__dict__.items():
if attr not in do_not_copy_by_ref:
new_group.__dict__[attr] = copy(value)
new_group._model = new
new.groups.append(new_group)
for group in self.groups:
new_group = new.groups.get_by_id(group.id)
# update awareness, as in the reaction copies
new_objects = []
for member in group.members:
if isinstance(member, Metabolite):
new_object = new.metabolites.get_by_id(member.id)
elif isinstance(member, Reaction):
new_object = new.reactions.get_by_id(member.id)
elif isinstance(member, Gene):
new_object = new.genes.get_by_id(member.id)
elif isinstance(member, Group):
new_object = new.groups.get_by_id(member.id)
else:
raise TypeError(
f"The group member {member!r} is unexpectedly not a "
f"metabolite, reaction, gene, nor another group."
)
new_objects.append(new_object)
new_group.add_members(new_objects)
try:
new._solver = deepcopy(self.solver)
# Cplex has an issue with deep copies
except Exception: # pragma: no cover
new._solver = copy(self.solver) # pragma: no cover
# it doesn't make sense to retain the context of a copied model so
# assign a new empty context
new._contexts = []
return new
[docs] def add_boundary(
self,
metabolite: Metabolite,
type: str = "exchange",
reaction_id: Optional[str] = None,
lb: Optional[float] = None,
ub: Optional[float] = None,
sbo_term: Optional[str] = None,
) -> Reaction:
"""
Add a boundary reaction for a given metabolite.
There are three different types of pre-defined boundary reactions:
exchange, demand, and sink reactions.
An exchange reaction is a reversible, unbalanced reaction that adds
to or removes an extracellular metabolite from the extracellular
compartment.
A demand reaction is an irreversible reaction that consumes an
intracellular metabolite.
A sink is similar to an exchange but specifically for intracellular
metabolites, i.e., a reversible reaction that adds or removes an
intracellular metabolite.
If you set the reaction `type` to something else, you must specify the
desired identifier of the created reaction along with its upper and
lower bound. The name will be given by the metabolite name and the
given `type`.
The change is reverted upon exit when using the model as a context.
Parameters
----------
metabolite : cobra.Metabolite
Any given metabolite. The compartment is not checked but you are
encouraged to stick to the definition of exchanges and sinks.
type : {"exchange", "demand", "sink"}
Using one of the pre-defined reaction types is easiest. If you
want to create your own kind of boundary reaction choose
any other string, e.g., 'my-boundary' (default "exchange").
reaction_id : str, optional
The ID of the resulting reaction. This takes precedence over the
auto-generated identifiers but beware that it might make boundary
reactions harder to identify afterwards when using `model.boundary`
or specifically `model.exchanges` etc. (default None).
lb : float, optional
The lower bound of the resulting reaction (default None).
ub : float, optional
The upper bound of the resulting reaction (default None).
sbo_term : str, optional
A correct SBO term is set for the available types. If a custom
type is chosen, a suitable SBO term should also be set (default None).
Returns
-------
cobra.Reaction
The created boundary reaction.
Examples
--------
>>> from cobra.io load_model
>>> model = load_model("textbook")
>>> demand = model.add_boundary(model.metabolites.atp_c, type="demand")
>>> demand.id
'DM_atp_c'
>>> demand.name
'ATP demand'
>>> demand.bounds
(0, 1000.0)
>>> demand.build_reaction_string()
'atp_c --> '
"""
ub = configuration.upper_bound if ub is None else ub
lb = configuration.lower_bound if lb is None else lb
types = {
"exchange": ("EX", lb, ub, sbo_terms["exchange"]),
"demand": ("DM", 0, ub, sbo_terms["demand"]),
"sink": ("SK", lb, ub, sbo_terms["sink"]),
}
if type == "exchange":
external = find_external_compartment(self)
if metabolite.compartment != external:
raise ValueError(
f"The metabolite is not an external metabolite (compartment is "
f"`{metabolite.compartment}` but should be `{external}`). "
f"Did you mean to add a demand or sink? If not, either change"
f" its compartment or rename the model compartments to fix this."
)
if type in types:
prefix, lb, ub, default_term = types[type]
if reaction_id is None:
reaction_id = f"{prefix}_{metabolite.id}"
if sbo_term is None:
sbo_term = default_term
if reaction_id is None:
raise ValueError(
"Custom types of boundary reactions require a custom "
"identifier. Please set the `reaction_id`."
)
if reaction_id in self.reactions:
raise ValueError(f"Boundary reaction '{reaction_id}' already exists.")
name = f"{metabolite.name} {type}"
rxn = Reaction(id=reaction_id, name=name, lower_bound=lb, upper_bound=ub)
rxn.add_metabolites({metabolite: -1})
if sbo_term:
rxn.annotation["sbo"] = sbo_term
self.add_reactions([rxn])
return rxn
[docs] def add_reactions(self, reaction_list: Iterable[Reaction]) -> None:
"""Add reactions to the model.
Reactions with identifiers identical to a reaction already in the
model are ignored.
The change is reverted upon exit when using the model as a context.
Parameters
----------
reaction_list : list
A list of `cobra.Reaction` objects
"""
def existing_filter(rxn: Reaction) -> bool:
"""Check if the reaction does not exists in the model.
Parameters
----------
rxn: cobra.Reaction
Returns
-------
bool
False if reaction exists, True if it doesn't.
If the reaction exists, will log a warning.
"""
if rxn.id in self.reactions:
logger.warning(f"Ignoring reaction '{rxn.id}' since it already exists.")
return False
return True
# First check whether the reactions exist in the model.
pruned = DictList(filter(existing_filter, reaction_list))
context = get_context(self)
# Add reactions. Also take care of genes and metabolites in the loop.
for reaction in pruned:
reaction._model = self
if context:
context(partial(setattr, reaction, "_model", None))
# Build a `list()` because the dict will be modified in the loop.
for metabolite in list(reaction.metabolites):
# TODO: Maybe this can happen with
# Reaction.add_metabolites(combine=False)
# TODO: Should we add a copy of the metabolite instead?
if metabolite not in self.metabolites:
self.add_metabolites(metabolite)
# A copy of the metabolite exists in the model, the reaction
# needs to point to the metabolite in the model.
else:
# FIXME: Modifying 'private' attributes is horrible.
stoichiometry = reaction._metabolites.pop(metabolite)
model_metabolite = self.metabolites.get_by_id(metabolite.id)
reaction._metabolites[model_metabolite] = stoichiometry
model_metabolite._reaction.add(reaction)
if context:
context(partial(model_metabolite._reaction.remove, reaction))
reaction.update_genes_from_gpr()
self.reactions += pruned
if context:
context(partial(self.reactions.__isub__, pruned))
# from cameo ...
self._populate_solver(pruned)
[docs] def remove_reactions(
self,
reactions: Union[str, Reaction, List[Union[str, Reaction]]],
remove_orphans: bool = False,
) -> None:
"""Remove reactions from the model.
The change is reverted upon exit when using the model as a context.
Parameters
----------
reactions : list or reaction or str
A list with reactions (`cobra.Reaction`), or their id's, to remove.
Reaction will be placed in a list. Str will be placed in a list and used to
find the reaction in the model.
remove_orphans : bool, optional
Remove orphaned genes and metabolites from the model as
well (default False).
"""
if isinstance(reactions, str) or hasattr(reactions, "id"):
warn("need to pass in a list")
reactions = [reactions]
context = get_context(self)
for reaction in reactions:
# Make sure the reaction is in the model
try:
reaction = self.reactions[self.reactions.index(reaction)]
except ValueError:
warn(f"{reaction} not in {self}")
else:
forward = reaction.forward_variable
reverse = reaction.reverse_variable
if context:
obj_coef = reaction.objective_coefficient
if obj_coef != 0:
context(
partial(
self.solver.objective.set_linear_coefficients,
{forward: obj_coef, reverse: -obj_coef},
)
)
context(partial(self._populate_solver, [reaction]))
context(partial(setattr, reaction, "_model", self))
context(partial(self.reactions.add, reaction))
self.remove_cons_vars([forward, reverse])
self.reactions.remove(reaction)
reaction._model = None
for met in reaction._metabolites:
if reaction in met._reaction:
met._reaction.remove(reaction)
if context:
context(partial(met._reaction.add, reaction))
if remove_orphans and len(met._reaction) == 0:
self.remove_metabolites(met)
for gene in reaction._genes:
if reaction in gene._reaction:
gene._reaction.remove(reaction)
if context:
context(partial(gene._reaction.add, reaction))
if remove_orphans and len(gene._reaction) == 0:
self.genes.remove(gene)
if context:
context(partial(self.genes.add, gene))
# remove reference to the reaction in all groups
associated_groups = self.get_associated_groups(reaction)
for group in associated_groups:
group.remove_members(reaction)
[docs] def add_groups(self, group_list: Union[str, Group, List[Group]]) -> None:
"""Add groups to the model.
Groups with identifiers identical to a group already in the model are
ignored.
If any group contains members that are not in the model, these members
are added to the model as well. Only metabolites, reactions, and genes
can have groups.
Parameters
----------
group_list : list or str or Group
A list of `cobra.Group` objects to add to the model. Can also be a single
group or a string representing group id. If the input is not a list, a
warning is raised.
"""
def existing_filter(new_group: Group) -> bool:
"""Check if the group does not exist.
Parameters
----------
new_group: cobra.Group
Group to check.
Returns
-------
bool
False if the group already exists, True if it doesn't.
"""
if new_group.id in self.groups:
logger.warning(
f"Ignoring group '{new_group.id}'" f" since it already exists."
)
return False
return True
if isinstance(group_list, str) or hasattr(group_list, "id"):
warn("need to pass in a list")
group_list = [group_list]
pruned = DictList(filter(existing_filter, group_list))
for group in pruned:
group._model = self
for member in group.members:
# If the member is not associated with the model, add it
if isinstance(member, Metabolite) and member not in self.metabolites:
self.add_metabolites([member])
if isinstance(member, Reaction) and member not in self.reactions:
self.add_reactions([member])
# TODO(midnighter): `add_genes` method does not exist.
# if isinstance(member, Gene):
# if member not in self.genes:
# self.add_genes([member])
self.groups += [group]
[docs] def remove_groups(self, group_list: Union[str, Group, List[Group]]) -> None:
"""Remove groups from the model.
Members of each group are not removed
from the model (i.e. metabolites, reactions, and genes in the group
stay in the model after any groups containing them are removed).
Parameters
----------
group_list : list or str or Group
A list of `cobra.Group` objects to remove from the model. Can also be a
single group or a string representing group id. If the input is not a list,
a warning is raised.
"""
if isinstance(group_list, str) or hasattr(group_list, "id"):
warn("need to pass in a list")
group_list = [group_list]
for group in group_list:
# make sure the group is in the model
if group.id not in self.groups:
logger.warning(f"{group!r} not in {self!r}. Ignored.")
else:
self.groups.remove(group)
group._model = None
[docs] def get_associated_groups(
self, element: Union[Reaction, Gene, Metabolite]
) -> List[Group]:
"""Get list of groups for element.
Returns a list of groups that an element (reaction, metabolite, gene)
is associated with.
Parameters
----------
element: `cobra.Reaction`, `cobra.Metabolite`, or `cobra.Gene`
Returns
-------
list of `cobra.Group`
All groups that the provided object is a member of
"""
# check whether the element is associated with the model
return [g for g in self.groups if element in g.members]
[docs] def add_cons_vars(
self, what: Union[List["CONS_VARS"], Tuple["CONS_VARS"]], **kwargs
) -> None:
"""Add constraints and variables to the model's mathematical problem.
Useful for variables and constraints that can not be expressed with
reactions and simple lower and upper bounds.
Additions are reversed upon exit if the model itself is used as
context.
Parameters
----------
what : list or tuple of optlang variables or constraints.
The variables or constraints to add to the model. Must be of
class `optlang.interface.Variable` or
`optlang.interface.Constraint`.
**kwargs : keyword arguments
Passed to solver.add()
"""
add_cons_vars_to_problem(self, what, **kwargs)
[docs] def remove_cons_vars(
self, what: Union[List["CONS_VARS"], Tuple["CONS_VARS"]]
) -> None:
"""Remove variables and constraints from problem.
Remove variables and constraints from the model's mathematical
problem.
Remove variables and constraints that were added directly to the
model's underlying mathematical problem. Removals are reversed
upon exit if the model itself is used as context.
Parameters
----------
what : list or tuple of optlang variables or constraints.
The variables or constraints to add to the model. Must be of
class `optlang.interface.Variable` or
`optlang.interface.Constraint`.
"""
remove_cons_vars_from_problem(self, what)
@property
[docs] def problem(self) -> "optlang.interface":
"""Get the interface to the model's underlying mathematical problem.
Solutions to cobra models are obtained by formulating a mathematical
problem and solving it. Cobrapy uses the optlang package to
accomplish that and with this property you can get access to the
problem interface directly.
Returns
-------
optlang.interface
The problem interface that defines methods for interacting with
the problem and associated solver directly.
"""
return self.solver.interface
@property
[docs] def variables(self) -> "Container":
"""Get the mathematical variables in the cobra model.
In a cobra model, most variables are reactions. However,
for specific use cases, it may also be useful to have other types of
variables. This property defines all variables currently associated
with the model's problem.
Returns
-------
optlang.container.Container
A container with all associated variables.
"""
return self.solver.variables
@property
[docs] def constraints(self) -> "Container":
"""Get the constraints in the cobra model.
In a cobra model, most constraints are metabolites and their
stoichiometries. However, for specific use cases, it may also be
useful to have other types of constraints. This property defines all
constraints currently associated with the model's problem.
Returns
-------
optlang.container.Container
A container with all associated constraints.
"""
return self.solver.constraints
@property
[docs] def boundary(self) -> List[Reaction]:
"""Boundary reactions in the model.
Reactions that either have no substrate or product.
Returns
-------
list
A list of reactions that either have no substrate or product and
only one metabolite overall.
"""
return [rxn for rxn in self.reactions if rxn.boundary]
@property
[docs] def exchanges(self) -> List[Reaction]:
"""Exchange reactions in model.
Reactions that exchange mass with the exterior. Uses annotations
and heuristics to exclude non-exchanges such as sink reactions.
Returns
-------
list
A list of reactions that satisfy the conditions for exchange reactions.
See Also
--------
cobra.medium.find_boundary_types
"""
return find_boundary_types(self, "exchange", None)
@property
[docs] def demands(self) -> List[Reaction]:
"""Demand reactions in model.
Irreversible reactions that accumulate or consume a metabolite in
the inside of the model.
Returns
-------
list
A list of reactions that are demand reactions (reactions that
accumulate/consume a metabolite irreversibly).
See Also
--------
cobra.medium.find_boundary_types
"""
return find_boundary_types(self, "demand", None)
@property
[docs] def sinks(self) -> List[Reaction]:
"""Sink reactions in model.
Reversible reactions that accumulate or consume a metabolite in
the inside of the model.
Returns
-------
list
A list of reactions that are demand reactions (reactions that
accumulate/consume a metabolite reversibly).
See Also
--------
cobra.medium.find_boundary_types
"""
return find_boundary_types(self, "sink", None)
[docs] def _populate_solver(
self,
reaction_list: List[Reaction],
metabolite_list: Optional[List[Metabolite]] = None,
) -> None:
"""Populate attached solver with constraints and variables.
Populate attached solver with constraints and variables that
model the provided reactions.
Parameters
----------
reaction_list: list
A list of cobra.Reaction to add to the solver. This list will be
constrained.
metabolite_list: list, optional
A list of cobra.Metabolite to add to the solver. This list will be
constrained (default None).
"""
constraint_terms = AutoVivification()
to_add = []
if metabolite_list is not None:
for met in metabolite_list:
to_add += [self.problem.Constraint(Zero, name=met.id, lb=0, ub=0)]
self.add_cons_vars(to_add)
for reaction in reaction_list:
if reaction.id not in self.variables:
forward_variable = self.problem.Variable(reaction.id)
reverse_variable = self.problem.Variable(reaction.reverse_id)
self.add_cons_vars([forward_variable, reverse_variable])
else:
reaction = self.reactions.get_by_id(reaction.id)
forward_variable = reaction.forward_variable
reverse_variable = reaction.reverse_variable
for metabolite, coeff in reaction.metabolites.items():
if metabolite.id in self.constraints:
constraint = self.constraints[metabolite.id]
else:
constraint = self.problem.Constraint(
Zero, name=metabolite.id, lb=0, ub=0
)
self.add_cons_vars(constraint, sloppy=True)
constraint_terms[constraint][forward_variable] = coeff
constraint_terms[constraint][reverse_variable] = -coeff
self.solver.update()
for reaction in reaction_list:
reaction = self.reactions.get_by_id(reaction.id)
reaction.update_variable_bounds()
for constraint, terms in constraint_terms.items():
constraint.set_linear_coefficients(terms)
[docs] def slim_optimize(
self, error_value: Optional[float] = float("nan"), message: Optional[str] = None
) -> float:
"""Optimize model without creating a solution object.
Creating a full solution object implies fetching shadow prices and
flux values for all reactions and metabolites from the solver
object. This necessarily takes some time and in cases where only one
or two values are of interest, it is recommended to instead use this
function which does not create a solution object returning only the
value of the objective. Note however that the `optimize()` function
uses efficient means to fetch values so if you need fluxes/shadow
prices for more than say 4 reactions/metabolites, then the total
speed increase of `slim_optimize` versus `optimize` is expected to
be small or even negative depending on how you fetch the values
after optimization.
Parameters
----------
error_value : float, None
The value to return if optimization failed due to e.g.
infeasibility. If None, raise `OptimizationError` if the
optimization fails (default float("nan")).
message : str, optional
Error message to use if the model optimization did not succeed (default
None).
Returns
-------
float
The objective value. Returns the error value if optimization failed and
error_value was not None.
Raises
------
OptimizationError
If error_value was set as None and the optimization fails.
"""
self.solver.optimize()
if self.solver.status == optlang.interface.OPTIMAL:
return self.solver.objective.value
elif error_value is not None:
return error_value
else:
assert_optimal(self, message)
[docs] def optimize(
self, objective_sense: Optional[str] = None, raise_error: bool = False
) -> "Solution":
"""
Optimize the model using flux balance analysis.
Parameters
----------
objective_sense : {None, 'maximize' 'minimize'}, optional
Whether fluxes should be maximized or minimized. In case of None,
the previous direction is used (default None).
raise_error : bool
If true, raise an OptimizationError if solver status is not
optimal (default False).
Returns
-------
Solution
Notes
-----
Only the most commonly used parameters are presented here. Additional
parameters for cobra.solvers may be available and specified with the
appropriate keyword argument.
"""
original_direction = self.objective.direction
self.objective.direction = {"maximize": "max", "minimize": "min"}.get(
objective_sense, original_direction
)
self.slim_optimize()
solution = get_solution(self, raise_error=raise_error)
self.objective.direction = original_direction
return solution
[docs] def repair(
self, rebuild_index: bool = True, rebuild_relationships: bool = True
) -> None:
"""Update all indexes and pointers in a model.
Parameters
----------
rebuild_index : bool
rebuild the indices kept in reactions, metabolites and genes
rebuild_relationships : bool
reset all associations between genes, metabolites, model and
then re-add them.
"""
if rebuild_index: # DictList indexes
self.reactions._generate_index()
self.metabolites._generate_index()
self.genes._generate_index()
self.groups._generate_index()
if rebuild_relationships:
for met in self.metabolites:
met._reaction.clear()
for gene in self.genes:
gene._reaction.clear()
for rxn in self.reactions:
rxn.update_genes_from_gpr()
for met in rxn._metabolites:
met._reaction.add(rxn)
# point _model to self
for dict_list in (self.reactions, self.genes, self.metabolites, self.groups):
for entity in dict_list:
entity._model = self
@property
[docs] def objective(self) -> Union[optlang.Objective]:
"""Get the solver objective.
With optlang, the objective is not limited to a simple linear summation of
individual reaction fluxes, making the return value ambiguous.
Henceforth, use `cobra.util.solver.linear_reaction_coefficients` to
get a dictionary of reactions with their linear coefficients (empty
if there are none).
"""
return self.solver.objective
@objective.setter
def objective(
self,
value: Union[
Dict[Reaction, "Container"],
str,
int,
Reaction,
"optlang.interface.Objective",
Basic,
],
) -> None:
"""Set the solver objective.
Parameters
----------
value: dict or str or int or Reaction or optlang.interface.Container, Reaction
or Basic
The set value can be dictionary (reactions as keys, linear coefficients as
values), string (reaction identifier), int (reaction index), Reaction or
problem.Objective or sympy expression directly interpreted as objectives.
When using in a context, this attribute can be set temporarily.
"""
if isinstance(value, Basic):
value = self.problem.Objective(value, sloppy=False)
if not isinstance(value, (dict, optlang.interface.Objective)):
try:
reactions = self.reactions.get_by_any(value)
except KeyError:
raise ValueError("invalid objective")
value = {rxn: 1 for rxn in reactions}
# TODO - check that it is reset with context.
set_objective(self, value, additive=False)
@property
[docs] def objective_direction(self) -> str:
"""Get the objective direction.
Returns
-------
str
Objective direction as string. Should be "max" or "min".
"""
return self.solver.objective.direction
@objective_direction.setter
@resettable
def objective_direction(self, value: str) -> None:
"""Set the objective direction.
When used in a context, this attribute is set temporarily.
Parameters
----------
value: {"max", "min"}
String of "max" or "min" for direction.
Raises
------
ValueError
If given direction isn't max or min.
"""
value = value.lower()
if value.startswith("max"):
self.solver.objective.direction = "max"
elif value.startswith("min"):
self.solver.objective.direction = "min"
else:
raise ValueError(f"Unknown objective direction '{value}'.")
[docs] def summary(
self,
solution: Optional["Solution"] = None,
fva: Union["pd.DataFrame", float, None] = None,
) -> "ModelSummary":
"""
Create a summary of the exchange fluxes of the model.
Parameters
----------
solution : cobra.Solution, optional
A previous model solution to use for generating the summary. If
``None``, the summary method will generate a parsimonious flux
distribution (default None).
fva : pd.DataFrame or float, optional
Whether or not to include flux variability analysis in the output.
If given, `fva` should either be a previous FVA solution matching the
model or a float between 0 and 1 representing the fraction of the
optimum objective to be searched (default None).
Returns
-------
cobra.ModelSummary
See Also
--------
Reaction.summary
Metabolite.summary
"""
from cobra.summary import ModelSummary
return ModelSummary(model=self, solution=solution, fva=fva)
[docs] def __enter__(self) -> "Model":
"""Record future changes to the model.
Record all future changes to the model, undoing them when a call to
__exit__ is received. Creates a new context and adds it to the stack.
Returns
-------
cobra.Model
Returns the model with context added.
"""
try:
self._contexts.append(HistoryManager())
except AttributeError:
self._contexts = [HistoryManager()]
return self
[docs] def __exit__(self, type, value, traceback) -> None:
"""Pop the top context manager and trigger the undo functions."""
context = self._contexts.pop()
context.reset()
[docs] def merge(
self,
right: "Model",
prefix_existing: Optional[str] = None,
inplace: bool = True,
objective: str = "left",
) -> "Model":
"""Merge two models to create a model with the reactions from both models.
Custom constraints and variables from right models are also copied
to left model, however note that, constraints and variables are
assumed to be the same if they have the same name.
Parameters
----------
right : cobra.Model
The model to add reactions from
prefix_existing : string or optional
Prefix the reaction identifier in the right that already exist
in the left model with this string (default None).
inplace : bool
Add reactions from right directly to left model object.
Otherwise, create a new model leaving the left model untouched.
When done within the model as context, changes to the models are
reverted upon exit (default True).
objective : {"left", "right", "sum"}
One of "left", "right" or "sum" for setting the objective of the
resulting model to that of the corresponding model or the sum of
both (default "left").
Returns
-------
cobra.Model
The merged model.
"""
if inplace:
new_model = self
else:
new_model = self.copy()
new_model.id = f"{self.id}_{right.id}"
new_reactions = deepcopy(right.reactions)
if prefix_existing is not None:
existing = new_reactions.query(lambda rxn: rxn.id in self.reactions)
for reaction in existing:
reaction.id = f"{prefix_existing}{reaction.id}"
new_model.add_reactions(new_reactions)
interface = new_model.problem
new_vars = [
interface.Variable.clone(v)
for v in right.variables
if v.name not in new_model.variables
]
new_model.add_cons_vars(new_vars)
new_cons = [
interface.Constraint.clone(c, model=new_model.solver)
for c in right.constraints
if c.name not in new_model.constraints
]
new_model.add_cons_vars(new_cons, sloppy=True)
new_model.objective = dict(
left=self.objective,
right=right.objective,
sum=self.objective.expression + right.objective.expression,
)[objective]
return new_model
[docs] def _repr_html_(self) -> str:
"""Get HTML represenation of the model.
Returns
-------
str
Model representation as HTML string.
"""
return f"""
<table>
<tr>
<td><strong>Name</strong></td>
<td>{self.id}</td>
</tr><tr>
<td><strong>Memory address</strong></td>
<td>{f"{id(self):x}"}</td>
</tr><tr>
<td><strong>Number of metabolites</strong></td>
<td>{len(self.metabolites)}</td>
</tr><tr>
<td><strong>Number of reactions</strong></td>
<td>{len(self.reactions)}</td>
</tr><tr>
<td><strong>Number of genes</strong></td>
<td>{len(self.genes)}</td>
</tr><tr>
<td><strong>Number of groups</strong></td>
<td>{len(self.groups)}</td>
</tr><tr>
<td><strong>Objective expression</strong></td>
<td>{format_long_string(str(self.objective.expression), 100)}</td>
</tr><tr>
<td><strong>Compartments</strong></td>
<td>{", ".join(v if v else k for k, v in
self.compartments.items())}</td>
</tr>
</table>"""