Exact Fastcore (MIP)
Example of implementation of the Fastcore algorithm with MIOM to extract a context specific model. Fastcore defines a set of core reactions that is forced to be active (carry a non-zero flux in steady state conditions), and minimizes the set of non-core reactions.
Cite
Vlassis, Pacheco, Sauter (2014). Fast reconstruction of compact context-specific metbolic network models. PLoS Comput. Biol. 10, e1003424.
The Fastcore algorithm is greedy approximation of the exact problem which can be modelled as a MIP problem. With MIOM, the exact problem can be defined in a few lines, using the opt_tol parameter to control the level of optimality required:
import miom
+ Exact Fastcore (MIP) - MIOM Exact Fastcore (MIP)
Example of implementation of the Fastcore algorithm with MIOM to extract a context specific model. Fastcore defines a set of core reactions that is forced to be active (carry a non-zero flux in steady state conditions), and minimizes the set of non-core reactions.
Cite
Vlassis, Pacheco, Sauter (2014). Fast reconstruction of compact context-specific metbolic network models. PLoS Comput. Biol. 10, e1003424.
The Fastcore algorithm is greedy approximation of the exact problem which can be modelled as a MIP problem. With MIOM, the exact problem can be defined in a few lines, using the opt_tol parameter to control the level of optimality required:
import miom
import numpy as np
# Use the flux-consistent subnetwork of the Human1 GEM model, available in the repository
@@ -30,7 +30,7 @@
Note that the algorithm requires that the metabolic network used is flux consistent. If blocked reactions are included in the core set, the MIP becomes infeasible, as those reactions cannot be selected with a non-zero flux in steady state conditions.
MIOM includes an implementation of the swiftcc algorithm to obtain flux consistent metabolic networks:
Cite
Tefagh, M., & Boyd, S. P. (2020). SWIFTCORE: a tool for the context-specific reconstruction of genome-scale metabolic networks. BMC bioinformatics, 21(1), 1-14. DOI: 10.1186/s12859-020-3440-y.
from miom.tools import consistent_subnetwork
# Get the consistent subnetwork and the lp problem solved by swiftcc
consistent, lp = consistent_subnetwork(m)
-
Flux Balance Analysis (LP)
A simple Flux Balance Analysis can be easily defined and solved with any of the commercial or open-source solvers available:
import miom
+ Flux Balance Analysis (LP) - MIOM Flux Balance Analysis (LP)
A simple Flux Balance Analysis can be easily defined and solved with any of the commercial or open-source solvers available:
import miom
# Load the iMM1865 mouse model from the repository
network = miom.mio.load_gem("@iMM1865")
target = "BIOMASS_reaction"
@@ -10,7 +10,7 @@
.get_fluxes(target))
print("Optimal flux is", flux)
-
iMAT: Integrative Metabolic Analysis Tool (MIP)
Example implementation of iMAT using MIOM. Note that this implementation supports also custom weights for the reactions. By default, weights for reactions in the "active" set are set to 1, and reactions in the "inactive" set are set to -1, so the objective function is exactly the same as the original iMAT.
Cite
Shlomi, T., Cabili, M. N., Herrgård, M. J., Palsson, B. Ø., & Ruppin, E. (2008). Network-based prediction of human tissue-specific metabolism. Nature biotechnology, 26(9), 1003-1010.
import miom
+ iMAT: Integrative Metabolic Analysis Tool (MIP) - MIOM iMAT: Integrative Metabolic Analysis Tool (MIP)
Example implementation of iMAT using MIOM. Note that this implementation supports also custom weights for the reactions. By default, weights for reactions in the "active" set are set to 1, and reactions in the "inactive" set are set to -1, so the objective function is exactly the same as the original iMAT.
Cite
Shlomi, T., Cabili, M. N., Herrgård, M. J., Palsson, B. Ø., & Ruppin, E. (2008). Network-based prediction of human tissue-specific metabolism. Nature biotechnology, 26(9), 1003-1010.
import miom
# Use the iHuman-GEM model
m = miom.load_gem('@SysBioChalmers/Human-GEM/v1.9.0')
@@ -23,7 +23,7 @@
)
.network)
print(m.num_reactions)
-
Sparse FBA (MIP)
The sparse FBA problem consists of finding the optimal flux value for a reaction minimizing the number of reactions with non-zero flux (l0-norm sparsity). Note that minimizing the l0-norm is a NP-hard problem, and so obtaing an optimal solution is not possible in many cases. The COBRA Toolbox includes different LP heuristics to minimize different approximations of the l0-norm. However, using a MIO formulation, it's is possible to obtain a solution which is close to optimal, with a lower number of reactions than the LP approximations, and very quickly.
Here is the implementation of sparse-FBA with MIOM:
import miom
+ Sparse FBA (MIP) - MIOM Sparse FBA (MIP)
The sparse FBA problem consists of finding the optimal flux value for a reaction minimizing the number of reactions with non-zero flux (l0-norm sparsity). Note that minimizing the l0-norm is a NP-hard problem, and so obtaing an optimal solution is not possible in many cases. The COBRA Toolbox includes different LP heuristics to minimize different approximations of the l0-norm. However, using a MIO formulation, it's is possible to obtain a solution which is close to optimal, with a lower number of reactions than the LP approximations, and very quickly.
Here is the implementation of sparse-FBA with MIOM:
import miom
# Load a genome-scale metabolic network. You can load SMBL or Matlab metabolic networks
# as well using the same method, but it requires to have the cobratoolbox python library
@@ -53,7 +53,7 @@
# Show the flux value of the biomass reaction
print("Optimal biomass flux:", V[m.get_reaction_id("BIOMASS_reaction")], "mmol/(h·gDW)")
-
Introduction
Warning
This library is functional but still in a very early stage. API is still not stable and might be changed in future versions.
MIOM (Mixed Integer Optimization for Metabolism) is a python library for creating and solving complex optimization problems using genome-scale metabolic networks, in just a few lines.
MIOM offers a high-level API that leverages the power of modern Mixed Integer Optimization (MIO) solvers to easily define steady-state metabolic optimization problems, from simple Flux Balance Analysis (FBA) simulations, to more complex problems, such as sparse FBA or context-specific reconstruction algorithms, and solve them the required level of optimality. It supports many free and commercial solvers thanks to the integration with the PICOS and the Python-MIP. It is also compatible and complementary to cobrapy.
Here is a simple example of the differences of implementing from scratch a simple Flux Balance Analysis (FBA), which is typically implemented as a Linear Program (LP) problem, and the Sparse FBA problem which is an Integer Programmming (IP) problem that solves the same FBA problem but minimizing the number of active reactions:
import miom
+ MIOM
Introduction
Warning
This library is functional but still in a very early stage. API is still not stable and might be changed in future versions.
MIOM (Mixed Integer Optimization for Metabolism) is a python library for creating and solving complex optimization problems using genome-scale metabolic networks, in just a few lines.
MIOM offers a high-level API that leverages the power of modern Mixed Integer Optimization (MIO) solvers to easily define steady-state metabolic optimization problems, from simple Flux Balance Analysis (FBA) simulations, to more complex problems, such as sparse FBA or context-specific reconstruction algorithms, and solve them the required level of optimality. It supports many free and commercial solvers thanks to the integration with the PICOS and the Python-MIP. It is also compatible and complementary to cobrapy.
Here is a simple example of the differences of implementing from scratch a simple Flux Balance Analysis (FBA), which is typically implemented as a Linear Program (LP) problem, and the Sparse FBA problem which is an Integer Programmming (IP) problem that solves the same FBA problem but minimizing the number of active reactions:
import miom
model = (miom
.load('@BiGG/Recon3D.miom')
.steady_state()
@@ -116,7 +116,7 @@
.network
)
print(fmc.num_reactions)
-
Advantages
- It's flexible: MIOM uses the PICOS and the Python-MIP libraries, which means you can use any solver supported by those libraries.
- It's easy to extend: MIOM is written in pure python, so you can easily extend it to solve more complex optimization problems.
- It makes the problem explicit: MIOM uses a declarative way to express the problem, so you can easily read and understand what you are solving and differences between algorithms.
- It's fast: MIOM leverages the power of MIO solvers to solve complex optimization problems. You can control the quality and speed of the solutions for your problem and get better solutions that the approximations (LP) of the original problem available in other constraint-based modeling libraries.
- It's lightweight: The library has a small number of dependencies, so it's easy to install and distribute also in HPC environments.
- It includes a compressed GEM format: MIOM can load and save the minimal information of the metabolic networks required for performing simulations into a compressed file compatible with numpy. The small size of the files allows you to quickly run online experiments so other people can reproduce your results. It also supports SBML and matlab formats if
cobratoolbox is installed. - It's open-source: MIOM is open-source and free to use. You can contribute to the development of MIOM by forking the repository and sending pull requests.
Documentation
The documentation of the library can be found at https://metexplore.github.io/miom/
How to cite
Manuscript in progress
License
GNU General Public License v3.0
mio.py
MiomNetwork
A minimal class to store a metabolic network.
Attributes:
Name Type Description Snumpy.ndarraystoichiometric matrix
Rnumpy.ndarrayStructured array with the reactions. The fields are: - id (int): reaction ID - lb (float): lower bound - ub (float): upper bound - subsystem (str): subsystem - gpr (str): gene-protein-reaction rule
Mnumpy.ndarrayStructured array with the metabolites. The fields are: - id (int): metabolite ID - name (str): metabolite name - formula (str): metabolite formula
Source code in miom/mio.py
class MiomNetwork:
- """A minimal class to store a metabolic network.
-
- Attributes:
- S (numpy.ndarray): stoichiometric matrix
- R (numpy.ndarray): Structured array with the reactions. The fields are:
- - id (int): reaction ID
- - lb (float): lower bound
- - ub (float): upper bound
- - subsystem (str): subsystem
- - gpr (str): gene-protein-reaction rule
- M (numpy.ndarray): Structured array with the metabolites. The fields are:
- - id (int): metabolite ID
- - name (str): metabolite name
- - formula (str): metabolite formula
- """
- def __init__(self, S, R, M):
- if S.shape[1] != R.shape[0]:
- raise ValueError("Num of reactions in the stoichiometric matrix "
- "are different from the number of reaction in R")
- self.S = S
- self.R = R
- self.M = M
-
- @property
- def num_reactions(self):
- """Number of reactions in the network.
-
- Returns:
- int: Number of reactions
- """
- return self.R.shape[0]
-
- @staticmethod
- def _find_reaction(rxn_id, R):
- if isinstance(rxn_id, Number):
- return rxn_id, R[rxn_id]
- for i, r in enumerate(R):
- if r['id'] == rxn_id:
- return i, r
- raise ValueError(f"Cannot find reaction {rxn_id}")
-
- def find_reaction(self, rxn_id):
- """Find a particular reaction in the metabolic network.
-
- Args:
- rxn_id (str): Name of the reaction
-
- Returns:
- numpy.ndarray: Structured array with the information of the reaction.
- """
- return MiomNetwork._find_reaction(rxn_id, self.R)
-
- def get_reaction_id(self, rxn):
- i, r = self.find_reaction(rxn)
- return i
-
- def find_reactions(self, rxn_ids):
- return [self.find_reaction(rxn_id) for rxn_id in rxn_ids]
-
- def find_reactions_from_pathway(self, pathway_name):
- return np.array([1 if pathway_name.lower() in rxn['subsystem'].lower() else 0 for rxn in self.R])
-
- def subnet(self, idxs):
- S_sub = self.S[:, idxs]
- R_sub = self.R[idxs]
- act_met = np.sum(np.abs(S_sub), axis=1) > 0
- M_sub = self.M[act_met]
- S_sub = S_sub[act_met, :]
- return MiomNetwork(S_sub, R_sub, M_sub)
-
- @property
- def object_size(self):
- bytes = self.S.__sizeof__() + self.R.__sizeof__() + self.M.__sizeof__()
- return bytes / 1024**2
-
num_reactions property readonly
Number of reactions in the network.
Returns:
Type Description intNumber of reactions
find_reaction(self, rxn_id)
Find a particular reaction in the metabolic network.
Parameters:
Name Type Description Default rxn_idstrName of the reaction
required
Returns:
Type Description numpy.ndarrayStructured array with the information of the reaction.
Source code in miom/mio.py
def find_reaction(self, rxn_id):
+ mio.py - MIOM mio.py
MiomNetwork
A minimal class to store a metabolic network.
Attributes:
Name Type Description Snumpy.ndarraystoichiometric matrix
Rnumpy.ndarrayStructured array with the reactions. The fields are: - id (int): reaction ID - lb (float): lower bound - ub (float): upper bound - subsystem (str): subsystem - gpr (str): gene-protein-reaction rule
Mnumpy.ndarrayStructured array with the metabolites. The fields are: - id (int): metabolite ID - name (str): metabolite name - formula (str): metabolite formula
num_reactions property readonly
Number of reactions in the network.
Returns:
Type Description intNumber of reactions
find_reaction(self, rxn_id)
Find a particular reaction in the metabolic network.
Parameters:
Name Type Description Default rxn_idstrName of the reaction
required
Returns:
Type Description numpy.ndarrayStructured array with the information of the reaction.
Source code in miom/mio.py
def find_reaction(self, rxn_id):
"""Find a particular reaction in the metabolic network.
Args:
@@ -100,7 +25,7 @@
# Store to file
with open(path_to_exported_file, 'wb') as f_out:
f_out.write(compressed)
-
load_gem(model_or_path)
Load a metabolic network from a file or URL.
The method supports any format supported by cobrapy (.xml, .yml, .json, .mat) or a miom compressed model (.miom) from a url or a local file path. For the cobra supported formats, you need the cobrapy package installed, and for .mat files, you need both cobrapy and scipy installed.
Parameters:
Name Type Description Default model_or_pathstrPath to a local file or URL pointing to the metabolic network. If the string starts with '@', the file will be loaded from the default github repository.
required
Returns:
Type Description MiomNetworkA MiomNetwork instance with the minimal information of the metabolic network required for simulations. It includes the stoichiometric matrix, the list of reactions with the lower and upper bounds, the associated genes and GPR rules, and the list of metabolites.
Source code in miom/mio.py
def load_gem(model_or_path):
+
load_gem(model_or_path)
Load a metabolic network from a file or URL.
The method supports any format supported by cobrapy (.xml, .yml, .json, .mat) or a miom compressed model (.miom) from a url or a local file path. For the cobra supported formats, you need the cobrapy package installed, and for .mat files, you need both cobrapy and scipy installed.
Parameters:
Name Type Description Default model_or_pathstrPath to a local file or URL pointing to the metabolic network. If the string starts with '@', the file will be loaded from the default github repository.
required
Returns:
Type Description MiomNetworkA MiomNetwork instance with the minimal information of the metabolic network required for simulations. It includes the stoichiometric matrix, the list of reactions with the lower and upper bounds, the associated genes and GPR rules, and the list of metabolites.
Source code in miom/mio.py
def load_gem(model_or_path):
"""Load a metabolic network from a file or URL.
The method supports any format supported by cobrapy (.xml, .yml, .json, .mat)
@@ -138,7 +63,7 @@
return cobra_to_miom(_read_cobra_model(file))
else:
return cobra_to_miom(model_or_path)
-
miom.py
BaseModel (ABC)
Base class for building LP/MIP optimization models using a high-level API for metabolic problems.
It implements the chainable methods to set-up a LP/MIP metabolic network problem. Two implementations are available: PicosModel and PythonMipModel. The PicosModel uses the PICOS library as a backend to interact with different solvers. The PythonMipModel uses the Python-MIP library to solve the model using the CBC or GUROBI solvers.
Note
Do not try to instantiate this class directly. Use the load() function instead. The method automatically selects the right implementation depending on the solver.
Source code in miom/miom.py
class BaseModel(ABC):
- """Base class for building LP/MIP optimization models
- using a high-level API for metabolic problems.
-
- It implements the chainable methods to set-up a LP/MIP
- metabolic network problem. Two implementations are available:
- PicosModel and PythonMipModel. The PicosModel uses the
- [PICOS library](https://picos-api.gitlab.io/picos/)
- as a backend to interact with different solvers. The PythonMipModel uses
- the [Python-MIP](https://www.python-mip.com/) library to solve the
- model using the CBC or GUROBI solvers.
-
- !!! note
- Do not try to instantiate this class directly. Use the [load()][miom.miom.load]
- function instead. The method automatically selects the right implementation
- depending on the solver.
-
- """
- def __init__(self, previous_step_model=None, miom_network=None, solver_name=None):
- self.previous_step_model = previous_step_model
- self.network = miom_network
- self.variables = None
- self.objective = None
- self._last_start_time = None
- self.last_solver_time = None
- if previous_step_model is not None:
- self._options = previous_step_model._options
- if miom_network is None:
- self.network = previous_step_model.network
- if solver_name is not None:
- self._options["solver"] = solver_name
- else:
- # Default recommended options
- self._options = {
- 'int_tol': 1e-8,
- 'feas_tol': 1e-8,
- 'opt_tol': 1e-5,
- 'verbosity': 0,
- 'solver': solver_name
- }
- self.problem = self.initialize_problem()
- self.setup(**self._options)
-
- @abstractmethod
- def initialize_problem(self):
- pass
-
- @abstractmethod
- def get_solver_status(self):
- pass
-
- @property
- def status(self):
- return self.get_solver_status()
-
-
- @_composable
- def setup(self, **kwargs):
- """Provide the options for the solver.
-
- Attributes:
- int_tol (float): Integrality tolerance for integer variables.
- Defaults to 1e-8.
- feas_tol (float): Feasibility tolerance. Defaults to 1e-8.
- opt_tol (float): Relative MIP gap tolerance for MIP problems.
- Defaults to 1e-5.
- verbosity (int): Values above 0 force the backends to be verbose.
- Use a value of 1 to show useful information about the search.
- Defaults to 0.
-
- Returns:
- BaseModel: the configured instance of the BaseModel
- """
- if self.problem is None:
- warnings.warn("Problem cannot be configured since it was not initialized")
- return False
- options = dict()
- int_tol = kwargs["int_tol"] if "int_tol" in kwargs else None
- feas_tol = kwargs["feas_tol"] if "feas_tol" in kwargs else None
- opt_tol = kwargs["opt_tol"] if "opt_tol" in kwargs else None
- verbosity = kwargs["verbosity"] if "verbosity" in kwargs else None
- solver = kwargs["solver"] if "solver" in kwargs else None
- if int_tol is not None:
- options["int_tol"] = int_tol
- if feas_tol is not None:
- options["feas_tol"] = feas_tol
- if opt_tol is not None:
- options["opt_tol"] = opt_tol
- if verbosity is not None:
- options["verbosity"] = verbosity
- if solver is not None:
- options["solver"] = solver
- self._options.update(options)
- # Calculate a min eps value to avoid numerical issues
- self._options["_min_eps"] = np.sqrt(2*self._options["int_tol"])
- return self._options
-
- @_composable
- def steady_state(self):
- """Add the required constraints for finding steady-state fluxes
-
- The method adds the $S * V = 0$ set of constraints, where $S$
- is the stoichiometric matrix and $V$ the flux variables.
-
- Returns:
- BaseModel: instance of BaseModel with the modifications applied.
- """
- pass
-
- @_composable
- def keep(self, reactions):
- """Force the inclusion of a list of reactions in the solution.
-
- Reactions have to be associated with positive weights in order to
- keep them in the final solution. Note that once keep() is called,
- the weights associated to the reactions selected to be kept will
- not be taken into account, as they will be forced to be kept in
- the solution.
-
- Args:
- reactions (list): List of reaction names, a binary vector
- indicating the reactions to keep, or a list of indexes
- with the reactions to keep.
-
- Returns:
- BaseModel: instance of BaseModel with the modifications applied.
- """
- if self.variables.indicators is None:
- raise ValueError("""No indicator variables for reactions,
- transform it to a subset selection problem calling
- subset_selection first, providing positive weights for the
- reactions you want to keep.""")
- if reactions is None or len(reactions) == 0:
- return False
- if isinstance(reactions, str):
- reactions = [reactions]
- # Check if it's a binary vector
- if len(reactions) == self.network.num_reactions and max(reactions)==1:
- reactions = np.where(reactions==1)[0]
- # Check if the reactions have an indicator variable
- reactions = set(self.network.find_reaction(rxn)[0] for rxn in reactions)
- available = set(rxn.index for rxn in self.variables.assigned_reactions if rxn.cost > 0)
- diff = reactions - available
- if len(diff) != 0:
- raise ValueError(f"""Only reactions associated with positive weights
- can be forced to be selected. The following reaction
- indexes have no indicator variables or are associated with
- negative weights: {diff}.""")
- valid_rxn_idx = reactions & available
- # Get the indexes of the indicator variables
- idxs = [i for i, r in enumerate(self.variables.assigned_reactions)
- if r.index in valid_rxn_idx]
- return dict(idxs=idxs, valid_rxn_idx=valid_rxn_idx)
-
- @_composable
- def subset_selection(self, rxn_weights, direction='max', eps=1e-2, strengthen=True):
- """Transform the current model into a subset selection problem.
-
- The subset selection problem consists of selecting the positive weighted
- reactions and remove the negative weighted reactions, subject to steady
- state constraints and/or additional constraints on fluxes, and maximizing
- the weighted sum of the (absolute) weights for the successfully selected reactions
- (with positive weights) and the successfully removed reactions (with negative
- weights). Selected reactions are forced to have an absolute flux value greater
- or equal to the threshold `eps` (1e-2 by default). Removed reactions should have a
- flux equal to 0.
-
- Each reaction is associated with a weight (positive or negative) provided
- in the parameter `rxn_weights`, and the objective is to select the reactions
- that optimizes the following expression:
-
- $$
- f(x) = \sum_i^n |w_i| * x_i
- $$
-
- where $x_i$ are the indicator variables for the reactions $i$ and $w_i$ are
- the weights for the reactions associated to the indicator variable. Indicator
- variables are automatically created for each reaction associated to a non-zero
- weight. Two (mutually exclusive) indicator variables are used for positive weighted
- reactions that are reversible to indicate whether there is positive or negative flux
- through the reaction. A single indicator variable is created for positive weighted
- non-reversible reactions, to indicate if the reaction is selected (has a non-zero
- flux greater or equal to `eps`) in which case the indicator variable is 1,
- or 0 otherwise.
-
- A single binary indicator variable is also created for negative weighted reactions,
- indicating whether the reaction was not selected (i.e, has 0 flux, in which case the
- indicator variable is 1) or not (in which case the indicator variable is 0).
-
- Args:
- rxn_weights (list): List of weights for each reaction. If a single value
- is provided, it is assumed to be the weight for all reactions.
- eps (float, optional): Min absolute flux value for weighted reactions
- to consider them active or inactive. Defaults to 1e-2.
-
- Returns:
- BaseModel: instance of BaseModel with the modifications applied.
- """
- # Calculate min valid EPS based on integrality tolerance
- min_eps = self._options["_min_eps"]
- if eps < min_eps:
- warnings.warn(f"Minimum epsilon value below min. allowed value, changed to {min_eps}.")
- eps = max(eps, min_eps)
- if not isinstance(rxn_weights, Iterable):
- rxn_weights = [rxn_weights] * self.network.num_reactions
- rxnw = _weighted_rxns(self.network.R, rxn_weights)
- if self.variables.indicators is None:
- self.variables._assigned_reactions = rxnw
- return dict(eps=eps, direction=direction, strengthen=strengthen)
- else:
- raise ValueError("The current model is already a subset selection problem.")
-
- @_composable
- def set_flux_bounds(self, rxn_id, min_flux=None, max_flux=None):
- """Change the flux bounds of a reaction.
-
- Args:
- rxn_id (str/int): name or id of the reaction to change
- min_flux (float, optional): Min flux value. Defaults to None.
- max_flux (float, optional): Max flux value. Defaults to None.
-
- Returns:
- BaseModel: instance of BaseModel with the modifications applied.
- """
- i, _ = self.network.find_reaction(rxn_id)
- return i
-
- @_composable
- def add_constraints(self, constraints):
- """Add a list of constraint to the model
-
- Args:
- constraints (list): List of expressions with the
- constraints.
-
- Returns:
- BaseModel: instance of BaseModel with the modifications applied.
- """
- for c in constraints:
- self.add_constraint(c)
- return len(constraints) > 0
-
- @_composable
- def exclude(self, indicator_values=None):
- """Exclude a solution from the set of feasible solutions.
-
- If the problem is a subset selection problem, it adds a new constraint
- to exclude the given values (or the current values of the indicator variables)
- from the set of feasible solutions. This is useful to force the solver to find
- alternative solutions in a manual fashion. https://doi.org/10.1371/journal.pcbi.1008730
-
-
- Args:
- indicator_values (list, optional): List of values for each indicator variable. Defaults to None.
-
- Returns:
- BaseModel: instance of BaseModel with the modifications applied.
- """
- if self.variables.indicators is None:
- raise ValueError("""The optimization model does not contain indicator variables.
- Make sure that the problem is a subset selection problem.""")
- if indicator_values is None:
- indicator_values = np.array(self.variables.indicator_values)
- return dict(indicator_values=indicator_values)
-
- @_composable
- def add_constraint(self, constraint):
- """Add a specific constraint to the model.
- The constraint should use existing variables already included in the model.
-
- Args:
- constraint: affine expression using model's variables.
- """
- pass
-
- @_composable
- def set_objective(self, cost_vector, variables, direction='max'):
- """Set the optmization objective of the model.
-
- Args:
- cost_vector (Iterable): List with the cost for each variable
- variables (Iterable): Variables used for the objective function
- direction (str, optional): Optimization direction (min or max). Defaults to 'max'.
- """
- if self.objective is not None:
- warnings.warn("Previous objective changed")
- self.objective = (cost_vector, variables, direction)
-
- def set_rxn_objective(self, rxn, direction='max'):
- """Set a flux objective
-
- Maximize or minimize the flux through a given reaction.
-
- Args:
- rxn (str): Name of the reaction to use for the optimization
- direction (str, optional): Minimization or maximization of
- the flux ('min' or 'max'). Defaults to 'max'.
-
- Returns:
- BaseModel: instance of BaseModel with the modifications applied.
- """
- i = self.network.get_reaction_id(rxn)
- cost = np.zeros(self.network.num_reactions)
- cost[i] = 1
- self.set_objective(cost, self.variables.fluxvars, direction=direction)
- return self
-
-
- def set_fluxes_for(self, reactions, tolerance=1e-6):
- warnings.warn("This method was renamed to fix_fluxes. It will dissappear in the v0.9.0", DeprecationWarning)
- return self.fix_fluxes(reactions, tolerance)
-
-
- def fix_fluxes(self, reactions, tolerance=1e-6):
- """Force the flux of certain reactions to match current values.
-
- After calling `.solve()` for a flux optimization problem (e.g. FBA), this
- method adds a new constraint to force the flux of the given reactions to
- match the current flux values found after the optimization.
-
- This is interesting for example to implement methods like sparse-FBA, where
- the optimization is no longer the flux but the number of active reactions,
- and a new constraint is required to preserve optimality of fluxes.
-
- Args:
- reactions (list): reaction or list of reactions
- tolerance (float, optional): Tolerance for the flux values
- (a solution is valid if the flux is within optimal value +/- tol.
- Defaults to 1e-6.
-
- Returns:
- BaseModel: instance of BaseModel with the modifications applied.
- """
- if isinstance(reactions, str):
- reactions = [reactions]
-
- for rid in reactions:
- if isinstance(rid, np.ndarray):
- # TODO: Test this path
- rid = rid['id']
- idx, rxn = self.network.find_reaction(rid)
- lb = max(rxn['lb'], self.variables.flux_values[idx] - tolerance)
- ub = min(rxn['ub'], self.variables.flux_values[idx] + tolerance)
- self.add_constraint(self.variables.fluxvars[idx] >= lb)
- self.add_constraint(self.variables.fluxvars[idx] <= ub)
- return self
-
- @_composable
- def reset(self):
- """Resets the original problem (removes all modifications)
-
- Returns:
- BaseModel: instance of BaseModel with the modifications applied.
- """
- if self.problem is None:
- warnings.warn("Problem is not initialized, nothing to reset")
- return False
- else:
- return True
-
- def obtain_subnetwork(
- self,
- extraction_mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
- comparator=Comparator.GREATER_OR_EQUAL,
- value=1e-8
- ):
- """Same as [select_subnetwork][miom.miom.BaseModel] but returns the network instead.
-
- See [select_subnetwork][miom.miom.BaseModel] for a detailed description of the method.
-
- Returns:
- MiomNetwork: A MiomNetwork with the selected subnetwork.
- """
- # If indicators are present and assigned,
- # take the subset of the network for which
- # the indicators of positive weighted reactions
- # are equal to 1
- if extraction_mode == ExtractionMode.ABSOLUTE_FLUX_VALUE:
- variables = self.variables.flux_values
- if variables is None:
- raise ValueError("The model does not contain flux variables. "
- "You need to call first steady_state() to add "
- "the flux variables")
- elif extraction_mode == ExtractionMode.INDICATOR_VALUE:
- variables = self.variables.indicator_values
- if variables is None:
- raise ValueError("The model does not contain indicator variables. "
- "You need to transform it to a subset selection problem "
- "by invoking subset_selection() first.")
- elif extraction_mode == ExtractionMode.REACTION_ACTIVITY:
- variables = self.variables.reaction_activity
- if variables is None:
- raise ValueError("The model does not contain reaction activity variables. "
- "You need to transform it to a subset selection problem "
- "by invoking subset_selection() first.")
- else:
- raise ValueError("Invalid extraction mode")
-
- if comparator == Comparator.GREATER_OR_EQUAL:
- selected = np.where(np.abs(variables) >= value)[0]
- elif comparator == Comparator.LESS_OR_EQUAL:
- selected = np.where(np.abs(variables) <= value)[0]
- else:
- raise ValueError("Invalid comparison")
- # TODO: Assigned reactions work only for indicator variables
- if extraction_mode == ExtractionMode.INDICATOR_VALUE:
- rxns = [self.variables.assigned_reactions[x] for x in selected]
- selected_idx = [rxn.index for rxn in rxns]
- elif extraction_mode == ExtractionMode.ABSOLUTE_FLUX_VALUE:
- selected_idx = selected
- else:
- raise NotImplementedError("Only indicator variables and absolute flux values are supported")
- return self.network.subnet(selected_idx)
-
- @_composable
- def select_subnetwork(
- self,
- mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
- comparator=Comparator.GREATER_OR_EQUAL,
- value=1e-8
- ):
- """Select a subnetwork and create a new BaseModel to operate on it.
-
- The new instance of the BaseModel is a new problem instance with no constraints.
- If the idea is to perform FBA simulations on this new subnetwork, remember to add
- the new constraints, especially the `steady_state`.
-
- Args:
- mode (ExtractionMode, optional): Method used to extract the subnetwork
- (based on flux values or using the indicator values).
- Defaults to ExtractionMode.ABSOLUTE_FLUX_VALUE.
- comparator (Comparator, optional): Comparator for the selected mode.
- Defaults to Comparator.GREATER_OR_EQUAL.
- value (float, optional): Value threshold for the mode and comparator selected.
- Defaults to 1e-8.
-
- Returns:
- [type]: [description]
- """
- return self.obtain_subnetwork(extraction_mode=mode,
- comparator=comparator,
- value=value)
-
- def get_values(self):
- """Get the values for the variables
-
- Returns:
- tuple: (V, X) where V are the flux values and X are the indicator values
- (if the problem is a MIP problem, for example if `subset_selection` was
- called)
- """
- return self.variables.values()
-
-
- def get_fluxes(self, reactions=None):
- """Get the flux values.
-
- Args:
- reactions (list, optional): Reaction or subset of reactions
- For which to obtain the flux values. Defaults to None.
-
- Returns:
- list: List with the flux values for all or the selected reactions.
- """
- if isinstance(reactions, str):
- return self.variables.flux_values[self.network.get_reaction_id(reactions)]
- if isinstance(reactions, Iterable):
- return {
- r['id']: self.variables.flux_values[self.network.get_reaction_id(r['id'])]
- for r in reactions
- }
- if reactions is None:
- return self.variables.flux_values
- else:
- raise ValueError("reactions should be an iterable of strings or a single string")
-
- @_composable
- def solve(self, verbosity=None, max_seconds=None):
- """Solve the current model and assign the values to the variables of the model.
-
- Args:
- verbosity (int, optional): Level of verbosity for the solver.
- Values above 0 will force the backend to show output information of the search. Defaults to None.
- max_seconds (int, optional): Max time in seconds for the search. Defaults to None.
- """
- self._last_start_time = perf_counter()
-
- @_composable
- def copy(self):
- pass
-
add_constraint(self, constraint)
Add a specific constraint to the model. The constraint should use existing variables already included in the model.
Parameters:
Name Type Description Default constraintaffine expression using model's variables.
required
Source code in miom/miom.py
@_composable
+ miom.py - MIOM miom.py
BaseModel
Base class for building LP/MIP optimization models using a high-level API for metabolic problems.
It implements the chainable methods to set-up a LP/MIP metabolic network problem. Two implementations are available: PicosModel and PythonMipModel. The PicosModel uses the PICOS library as a backend to interact with different solvers. The PythonMipModel uses the Python-MIP library to solve the model using the CBC or GUROBI solvers.
Note
Do not try to instantiate this class directly. Use the load() function instead. The method automatically selects the right implementation depending on the solver.
add_constraint(self, constraint)
Add a specific constraint to the model. The constraint should use existing variables already included in the model.
Parameters:
Name Type Description Default constraintaffine expression using model's variables.
required
Source code in miom/miom.py
@_composable
def add_constraint(self, constraint):
"""Add a specific constraint to the model.
The constraint should use existing variables already included in the model.
@@ -640,7 +150,7 @@
idxs = [i for i, r in enumerate(self.variables.assigned_reactions)
if r.index in valid_rxn_idx]
return dict(idxs=idxs, valid_rxn_idx=valid_rxn_idx)
-
obtain_subnetwork(self, extraction_mode=<ExtractionMode.ABSOLUTE_FLUX_VALUE: 'flux_value'>, comparator=<Comparator.GREATER_OR_EQUAL: 'geq'>, value=1e-08)
Same as select_subnetwork but returns the network instead.
See select_subnetwork for a detailed description of the method.
Returns:
Type Description MiomNetworkA MiomNetwork with the selected subnetwork.
Source code in miom/miom.py
def obtain_subnetwork(
+
obtain_subnetwork(self, extraction_mode=<ExtractionMode.ABSOLUTE_FLUX_VALUE: 'flux_value'>, comparator=<Comparator.GREATER_OR_EQUAL: 'geq'>, value=1e-08)
Same as select_subnetwork but returns the network instead.
See select_subnetwork for a detailed description of the method.
Returns:
Type Description MiomNetworkA MiomNetwork with the selected subnetwork.
Source code in miom/miom.py
def obtain_subnetwork(
self,
extraction_mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
comparator=Comparator.GREATER_OR_EQUAL,
@@ -895,18 +405,7 @@
return dict(eps=eps, direction=direction, strengthen=strengthen)
else:
raise ValueError("The current model is already a subset selection problem.")
-
Comparator (str, Enum)
Comparator enum to use with select_subnetwork()
Attributes:
Name Type Description GREATER_OR_EQUALstrSelect those variables with a value greater or equal than the one provided.
LESS_OR_EQUALstrSelect those variables with a value less or equal than the one provided.
Source code in miom/miom.py
class Comparator(str, Enum):
- """Comparator enum to use with [select_subnetwork()][miom.miom.BaseModel]
-
- Attributes:
- GREATER_OR_EQUAL (str): Select those variables
- with a value greater or equal than the one provided.
- LESS_OR_EQUAL (str): Select those variables
- with a value less or equal than the one provided.
- """
- GREATER_OR_EQUAL = "geq",
- LESS_OR_EQUAL = "leq"
-
ExtractionMode (str, Enum)
Method to select the subnetwork to be extracted.
This mode works with the method select_subnetwork only after a subset selection problem was succesfully solved. If the model is configured as a subset selection problem, the user can extract a subnetwork after the method solve() was called. The selection of the subnetwork can be done using the indicator variables or the flux variables.
For more information, please read the documentation of the method subset_selection
Attributes:
Name Type Description ABSOLUTE_FLUX_VALUEstrUse a decision criterion based on the value of the fluxes. For example, by selecting all the reactions that have an absolute flux value above a certain threshold.
INDICATOR_VALUEstrUse the binary indicator variables to select the subnetwork. Binary indicator variables are created for a subset selection problem (after calling subset_selection). Two indicators are created for each positive weighted and reversible reaction, to indicate if there is a non-zero positive flux or negative flux respectively. A single binary indicator variable is created for each negative weighted reaction. In this case, an indicator value of 1 indicates that the reaction was succesfully removed, and 0 otherwise. You can use the value of the indicator variables to select a subnetwork after solving. For example, if all the reactions have a negative weight (since the goal is, for example, to minimize the number of reactions subject to some other constraints) and you want to select only the reactions that were not removed after solving, you can use the indicator variables with a value of 0.
Usage example:
m.
+
Comparator
Comparator enum to use with select_subnetwork()
Attributes:
Name Type Description GREATER_OR_EQUALstrSelect those variables with a value greater or equal than the one provided.
LESS_OR_EQUALstrSelect those variables with a value less or equal than the one provided.
ExtractionMode
Method to select the subnetwork to be extracted.
This mode works with the method select_subnetwork only after a subset selection problem was succesfully solved. If the model is configured as a subset selection problem, the user can extract a subnetwork after the method solve() was called. The selection of the subnetwork can be done using the indicator variables or the flux variables.
For more information, please read the documentation of the method subset_selection
Attributes:
Name Type Description ABSOLUTE_FLUX_VALUEstrUse a decision criterion based on the value of the fluxes. For example, by selecting all the reactions that have an absolute flux value above a certain threshold.
INDICATOR_VALUEstrUse the binary indicator variables to select the subnetwork. Binary indicator variables are created for a subset selection problem (after calling subset_selection). Two indicators are created for each positive weighted and reversible reaction, to indicate if there is a non-zero positive flux or negative flux respectively. A single binary indicator variable is created for each negative weighted reaction. In this case, an indicator value of 1 indicates that the reaction was succesfully removed, and 0 otherwise. You can use the value of the indicator variables to select a subnetwork after solving. For example, if all the reactions have a negative weight (since the goal is, for example, to minimize the number of reactions subject to some other constraints) and you want to select only the reactions that were not removed after solving, you can use the indicator variables with a value of 0.
Usage example:
m.
steady_state().
add_constraints(...).
# Convert to a subset selection, where each reaction
@@ -928,123 +427,10 @@
comparator=Comparator.LESS_OR_EQUAL,
value=0.5
).network
-
Source code in miom/miom.py
class ExtractionMode(str, Enum):
- """Method to select the subnetwork to be extracted.
-
- This mode works with the method [select_subnetwork][miom.miom.BaseModel]
- only after a subset selection problem was succesfully solved.
- If the model is configured as a subset selection problem, the user can
- extract a subnetwork after the method solve() was called. The selection
- of the subnetwork can be done using the indicator variables or the flux
- variables.
-
- For more information, please read the documentation of the method
- [subset_selection][miom.miom.BaseModel]
-
- Attributes:
- ABSOLUTE_FLUX_VALUE (str): Use a decision criterion based on the
- value of the fluxes. For example, by selecting all the reactions
- that have an absolute flux value above a certain threshold.
-
- INDICATOR_VALUE (str): Use the binary indicator variables to select
- the subnetwork. Binary indicator variables are created for a
- subset selection problem (after calling [subset_selection][miom.miom.BaseModel]).
- Two indicators are created for each positive weighted and reversible reaction,
- to indicate if there is a non-zero positive flux or negative flux respectively.
- A single binary indicator variable is created for each negative weighted reaction.
- In this case, an indicator value of 1 indicates that the reaction was succesfully
- removed, and 0 otherwise. You can use the value of the indicator variables to
- select a subnetwork after solving. For example, if all the reactions have a
- negative weight (since the goal is, for example, to minimize the number of
- reactions subject to some other constraints) and you want to select only the
- reactions that were not removed after solving, you can use the indicator
- variables with a value of 0.
-
- Usage example:
- ```python
- m.
- steady_state().
- add_constraints(...).
- # Convert to a subset selection, where each reaction
- # has a weight of -1.
- subset_selection(-1).
- # Solve the optimization problem
- .solve()
- # Get the subnetwork selecting the reactions
- # with an indicator value below 0.5. Since variables
- # are binary, it basically selects all the reactions
- # associated with a binary indicator value of 0.
- # This corresponds to the reactions that were not
- # removed after solving the problem (since all the
- # reactions have negative weights, and their binary
- # indicator variables are 1 if they were successfully
- # removed, and 0 if they could not be removed).
- .select_subnetwork(
- mode=ExtractionMode.INDICATOR_VALUE,
- comparator=Comparator.LESS_OR_EQUAL,
- value=0.5
- ).network
- ```
- """
- ABSOLUTE_FLUX_VALUE = "flux_value",
- INDICATOR_VALUE = "indicator_value",
- REACTION_ACTIVITY = "reaction_activity"
-
Solvers (str, Enum)
Solvers supported by the miom module.
Please refer to https://picos-api.gitlab.io/picos/introduction.html to see the list of supported solvers using the PICOS backend. The Python-MIP backend only supports the GUROBI and CBC solvers.
Note that in some cases, the Python-MIP backend (GUROBI/CBC) might be faster setting up the problem than the PICOS backend since it has less overhead.
Attributes:
Name Type Description GUROBI_PYMIPstrRecommended solver for most problems (LP, MIP). It uses the Python-MIP backend. Note that GUROBI is a commercial software which require a license. Free academic licenses are also available at https://gurobi.com/free/.
GUROBIstrFor using GUROBI with the PICOS backend instead of Python-MIP.
COIN_OR_CBCstrFree LP/MIP solver with good performance, provided with Python-MIP.
CPLEXstrCommercial solver with a performance similar to GUROBI. Only supported by the PICOS backend.
SCIPstrFree academic solver, supported by the PICOS backend.
GLPKstrOpen source solver, supported by the PICOS backend.
MOSEKstrCommercial solver, supported by the PICOS backend.
Source code in miom/miom.py
class Solvers(str, Enum):
- """Solvers supported by the miom module.
-
- Please refer to https://picos-api.gitlab.io/picos/introduction.html to see
- the list of supported solvers using the PICOS backend. The Python-MIP backend
- only supports the GUROBI and CBC solvers.
-
- Note that in some cases, the Python-MIP backend (GUROBI/CBC) might be faster
- setting up the problem than the PICOS backend since it has less overhead.
-
- Attributes:
- GUROBI_PYMIP (str): Recommended solver for most problems (LP, MIP).
- It uses the Python-MIP backend. Note that GUROBI is a commercial
- software which require a license. Free academic licenses are
- also available at https://gurobi.com/free/.
- GUROBI (str): For using GUROBI with the PICOS backend instead of Python-MIP.
- COIN_OR_CBC (str): Free LP/MIP solver with good performance, provided with Python-MIP.
- CPLEX (str): Commercial solver with a performance similar to GUROBI. Only
- supported by the PICOS backend.
- SCIP (str): Free academic solver, supported by the PICOS backend.
- GLPK (str): Open source solver, supported by the PICOS backend.
- MOSEK (str): Commercial solver, supported by the PICOS backend.
- """
- GUROBI_PYMIP = "gurobi_pymip",
- GUROBI = "gurobi",
- COIN_OR_CBC = "cbc",
- CPLEX = "cplex",
- GLPK = "glpk",
- SCIP = "scip",
- CVXOPT = "cvxopt",
- MOSEK = "mosek"
-
Status (str, Enum)
Subset of PICOS solver status codes.
This subset of states offers compatibility between Python-MIP and PICOS. Status codes can be obtained after solving an optimization problem with the property status. For example:
>>> import miom
+
Solvers
Solvers supported by the miom module.
Please refer to https://picos-api.gitlab.io/picos/introduction.html to see the list of supported solvers using the PICOS backend. The Python-MIP backend only supports the GUROBI and CBC solvers.
Note that in some cases, the Python-MIP backend (GUROBI/CBC) might be faster setting up the problem than the PICOS backend since it has less overhead.
Attributes:
Name Type Description GUROBI_PYMIPstrRecommended solver for most problems (LP, MIP). It uses the Python-MIP backend. Note that GUROBI is a commercial software which require a license. Free academic licenses are also available at https://gurobi.com/free/.
GUROBIstrFor using GUROBI with the PICOS backend instead of Python-MIP.
COIN_OR_CBCstrFree LP/MIP solver with good performance, provided with Python-MIP.
CPLEXstrCommercial solver with a performance similar to GUROBI. Only supported by the PICOS backend.
SCIPstrFree academic solver, supported by the PICOS backend.
GLPKstrOpen source solver, supported by the PICOS backend.
MOSEKstrCommercial solver, supported by the PICOS backend.
Status
Subset of PICOS solver status codes.
This subset of states offers compatibility between Python-MIP and PICOS. Status codes can be obtained after solving an optimization problem with the property status. For example:
>>> import miom
>>> miom.load('@iMM1865').steady_state().set_rxn_objective('BIOMASS_reaction').solve().status
{'status': 'optimal', 'objective_value': 798.811, 'elapsed_seconds': 0.479}
-
Source code in miom/miom.py
class Status(str, Enum):
- """Subset of PICOS solver status codes.
-
- This subset of states offers compatibility between Python-MIP and PICOS.
- Status codes can be obtained after solving an optimization problem with
- the property `status`. For example:
-
- ```python
- >>> import miom
- >>> miom.load('@iMM1865').steady_state().set_rxn_objective('BIOMASS_reaction').solve().status
- {'status': 'optimal', 'objective_value': 798.811, 'elapsed_seconds': 0.479}
- ```
- """
- OPTIMAL = "optimal",
- FEASIBLE = "feasible"
- INFEASIBLE = "infeasible",
- UNBOUNDED = "unbounded",
- ILLPOSED = "illposed",
- EMPTY = "empty",
- PREMATURE = "premature",
- UNKNOWN = "unknown"
-
load(network, solver=None)
Create a MIOM optimization model for a given solver. If the solver is Coin-OR CBC, an instance of PythonMipModel is used (which uses the Python-MIP lib as the backend). Otherwise, a PicosModel is created (which uses PICOS as the backend).
Examples:
Example of how to perform FBA to maximize flux through the BIOMASS_reaction in the iMM1865 model:
>>> import miom
+
load(network, solver=None)
Create a MIOM optimization model for a given solver. If the solver is Coin-OR CBC, an instance of PythonMipModel is used (which uses the Python-MIP lib as the backend). Otherwise, a PicosModel is created (which uses PICOS as the backend).
Examples:
Example of how to perform FBA to maximize flux through the BIOMASS_reaction in the iMM1865 model:
>>> import miom
>>> network = miom.load_gem("@iMM1865")
>>> flx = (miom
.load(network)
@@ -1052,7 +438,7 @@
.set_rxn_objective("BIOMASS_reaction")
.solve(verbosity=1)
.get_fluxes('BIOMASS_reaction'))
-
Parameters:
Name Type Description Default networkmiom_network or strA miom metabolic network or a valid file that can be imported with load_gem().
required solverSolverThe solver to be used. Defaults to Solver.GLPK.
None
Returns:
Type Description BaseModelA BaseModel object, which can be PythonMipModel if CBC solver is used, or a PicosModel otherwise.
Source code in miom/miom.py
def load(network, solver=None):
+
Parameters:
Name Type Description Default networkmiom_network or strA miom metabolic network or a valid file that can be imported with load_gem().
required solverSolverThe solver to be used. Defaults to Solver.GLPK.
None
Returns:
Type Description BaseModelA BaseModel object, which can be PythonMipModel if CBC solver is used, or a PicosModel otherwise.
Source code in miom/miom.py
def load(network, solver=None):
"""
Create a MIOM optimization model for a given solver.
If the solver is Coin-OR CBC, an instance of PythonMipModel is used
@@ -1105,7 +491,7 @@
raise Exception("""PICOS is not installed. Please install it with pip,
or reinstall miom with the right options, for example with:
'pip install miom[all]'.""")
-