Model optimization

Optimize

Base class for model optimizations. Supporting COBRApy and gurobipy interfaces.

class f2xba.optimization.optimize.Optimize(model_type, fname, cobra_model=None)[source]

Support optimization of SBML encoded extended models generated by f2xba.

Optimize is parent to FbaOptimization, EcmOptimization and RbaOptimization classes. Optimize provides common functionality, including support for thermodynamics constraint model variants.

Both COBRApy and gurobipy interfaces are supported. Use of gurobipy interface requires GUROBI and gurobipy being installed on your system. The gurobipy interface has significant speed advantages, especially for larger amd more complex models.

Example 1: Load and optimize GECKO model using the COBRApy interface.

import cobra
from f2xba import EcmOptimization

ecm = cobra.io.read_sbml_model('iML1515_GECKO.xml')
eo = EcmOptimization('iML1515_GECKO.xml', ecm)

ecm.medium = {rid: 1000.0 for rid in medium}
solution = ecm.optimize()

Example 2: Load and optimize GECKO model using the gurobipy interface.

from f2xba import EcmOptimization

eo = EcmOptimization('iML1515_GECKO.xml')
eo.medium = {rid: 1000.0 for rid in medium}
solution = eo.optimize()
m_dict

Information extracted from the SBML encoded model.

avg_enz_saturation

Average enzyme saturation level configured in ECM and RBA models.

protein_per_gdw

Total protein per gram cell dry weight (gP/gDW) configured in ECM models.

modeled_protein_mf

Protein mass fraction explicitly modelled in the ECM.

importer_km

Default Michaelis Menten constant KM (mmol/l) configured in the RBA model.

gpm

Instance of gurobipy.Model, in case gurobipy interface is used.

locus2uid

Map gene identifier to UniProt protein identifier.

id2groups

Map reaction identifier to group names configured in the SBML model.

rdata_model

Information related to reactions.

uptake_rids

Exchange reactions used for exchange of metabolites with the environment.

rdata

Information related to reactions (prefixes R_ and M_ removed).

net_rdata

Information related to net reactions (iso-reactions and fwd/rev directions combined).

rids_catalyzed

Map reaction identifier to related genes identifiers.

report_model_size()[source]

Print information about the model size.

eo.report_model_size()
get_tx_metab_genes()[source]

Get gene identifiers grouped by transport and metabolic reactions.

Genes required in both transport and metabolic enzymes are designated transport related genes.

Example: Extract sets of transport and metabolic related genes.

tx_genes, metab_genes = eo.get_tx_metab_genes()
Returns:

gene identifiers grouped by transport and metabolic enzymes

Return type:

(set(str), set(str))

get_rids_blocked_by(genes)[source]

Identify reactions that get blocked when deleting selected genes.

Example: Identify reactions that are blocked by the strain background.

strain_background = {'b0063', 'b0062', 'b0061', 'b3904', 'b3903', 'b3902', 'b0344'}
blocked_rids = eo.get_rids_blocked_by(strain_background)
Parameters:

genes (str or list(str)) – gene or genes

Returns:

reaction ids affected by gene deletions

Returns:

list(str)

get_variable_bounds(var_ids)[source]

Retrieve upper/lower bounds of selected model variables.

Bounds can be queried for a single variable or a list of variables.

Example: Retrieve variable bounds for single variable, list of variables.

bounds_1 = eo.get_variable_bounds('V_PC_total')
bounds_2 = eo.get_variable_bounds(['POR5_iso1', 'POR5_iso1_REV'])
Parameters:

var_ids (str or list(str)) – variable identifier or a list of variable identifiers

Returns:

variable ids with respective lower und upper bounds

Return type:

dict(str, tuple)

set_variable_bounds(variable_bounds)[source]

Set upper/lower bounds of selected model variables.

Bound values set to ‘None’ will not be modified.

Example: Block iso-reaction POR5_iso1 in forward and reverse directions. Only upper bounds need to be changed to 0.0.

orig_bounds = eo.set_variable_bounds({'POR5_iso1': (None, 0.0), 'POR5_iso1_REV': (None, 0.0)})
Parameters:

variable_bounds (dict(str, tuple)) – variable identifiers with tuple of lower and upper bounds

Returns:

original bounds

Return type:

dict(str, tuple)

get_objective()[source]

Retrieve model optimization objective.

Example: Query optimization objective of the model.

eo.get_objective()
Returns:

variable coefficients and optimization direction

Return type:

dict(str, float), str

set_objective(objective, direction='max')[source]

Set optimization objective.

Example: Modify the optimization objective of a model.

eo.set_objective({'BIOMASS_Ec_iML1515_WT_75p37M': 1.0}, 'max')
Parameters:

  • objective (dict(str, float)) – variable identifiers with coefficients

  • direction (str) – direction of optimization: ‘min’ or ‘max’ (default: ‘max’)

set_tfa_metab_concentrations(metab_concs)[source]

Configure metabolite concentration ranges in mol/l for TD constraint models.

Metabolomics data for a given condition can be used to further constrain models that contain thermodynamics constraints. Such models have optimization variables for species concentrations. Minimum and maximum concentrations for selected species can be configured using this method. Concentration bounds will not be set for ‘None’ values.

Example: Set minimum and maximum concentration values (mol/l) for selected metabolites.

molar_concs = {'atp_c': (0.00178, 0.00712), 'adp_c': (5.8e-05, 0.00023), 'datp_c': (9.05e-05, 0.00036)}
eo.set_tfa_metab_concentrations(molar_concs)
Parameters:

metab_concs (dict(str, tuple)) – min and max concentrations for species in mol/l

Returns:

original concentrations (min, max)

Return type:

dict (str, tuple(float))

update_relaxation(fluxes, eps=0.5)[source]

TFA slack model relaxation for a specified flux distribution.

In case optimizations of a model with added thermodynamics constraints leads to infeasable solutions, model relaxation is required. Model relaxation is performed on a TFA model with slack variables configured for ∆rG’˚ (standard transformed Gibbs energy of reaction) variables. Use TfaModel.export_slack_model() to generate such a model. This method is used to determine adjusted ∆rG’˚ variable bounds that make the TD model feasible and to configure the updated bounds to support optimization of further conditions.

Example 1: Perform TFA model relaxation using COBRApy interface. With the additional constraint of achieving at least 95% of the growth rate of the model without TD constraints.

import cobra
from f2xba import FbaOptimization

tfa_slack_model = cobra.io.read_sbml_model('iML1515_TFA_slack.xml')
fo = FbaOptimization('iML1515_TFA_slack.xml', tfa_slack_model)

tfa_slack_model.reactions.get_by_id('BIOMASS_Ec_iML1515_core_75p37M').lower_bound = 0.95 * wt_gr
tfa_slack_model.medium = medium

solution = tfa_slack_model.optimize()
drg0_relaxations = fo.update_relaxation(solution.fluxes)

Example 2: Perform TFA model relaxation using cobrapy interface.

from f2xba import FbaOptimization

fo = FbaOptimization('iML1515_TFA_slack.xml')

orig_bounds = fo.set_variable_bounds({'BIOMASS_Ec_iML1515_core_75p37M':(0.95 * wt_gr, None)})
fo.medium = medium

solution = fo.optimize()
fo.set_variable_bounds(orig_bounds)
drg0_relaxations = fo.update_relaxation(solution.fluxes)
Parameters:

  • fluxes (dict) – fluxes of a TFA slack model optimization solution

  • eps (float) – additional margin in kJ/mol (default: 0.5)

Returns:

∆rG’˚ variables with updated bounds.

Return type:

dict(str, dict(str,float)

modify_stoic(constr_id, var_id, new_val)[source]

Modify a stoichiometric coefficient using gurobipy interface.

Stoichiometric coefficients in a COBRApy model can be modified using the COBRApy method cobra.core.reaction.add_metabolites(). This method supports the modification of a stoichiometric coefficient when using the gurobipy interface. As model contexts are not supported, the coefficient should be reset to the original value after optimization.

Example: Remove the cofactor requirement of hemeA in anaerobic conditions.

if 'EX_o2_e' not in medium:
    orig_stoic = eo.modify_stoic('M_hemeA_c', 'R_CofactorSynth', 0.0)

eo.optimize()

if 'EX_o2_e' not in medium:
    eo.modify_stoic('M_hemeA_c', 'R_CofactorSynth', orig_stoic)
Parameters:

  • constr_id (str) – constraint identifier, e.g. species id

  • var_id (str) – variable identifier, e.g. reaction id

  • new_val (float) – stoichiometric coefficient (negative for consumation)

Returns:

original stoichiometric coefficient

Return type:

float

optimize(alt_model=None)[source]

Optimize gurobipy model and return COBRApy compliant solution.

Example: Optimize the model.

eo.optimize()

Thermodynamics constrained models are supported inline with pyTFA package.

Parameters:

alt_model (gurobipy.Model) – gurobipy model (default: None)

Returns:

optimization solution

Return type:

Solution

pfba(fraction_of_optimum=1.0)[source]

Perform parsimonious FBA optimization (pFBA) on a FBA model using the gurobipy interface.

Example: Execute a pFBA optimization using the gurobipy interface.

solution = eo.pfba()
Parameters:

fraction_of_optimum (float) – factor to scale FBA objective (default: 1.0)

Returns:

optimization solution

Returns:

Solution

fva(rids=None, fraction_of_optimum=1.0)[source]

Perform a Flux Variability Analysis (FVA) on a FBA model using the gurobipy interface.

A list of reaction ids can be provided to limit the scope of the FVA analysis.

Example: Execute FVA analyis for a set of reactions using the gurobipy interface.

eo.fva(['CS', 'ME1', 'POR5_iso1', 'POR5_iso1_REV'])
Parameters:

  • rids (list(str)) – reaction identifiers (without leading R_) (default: None)

  • fraction_of_optimum (float) – scaling of wild type objective value (default: 1.0)

Returns:

table with minimum and maximum reaction fluxes

Return type:

pandas.DataFrame

gene_knock_outs(genes)[source]

Block reactions that are catalyzed by specified gene/genes.

Example: Block reactions that depend on the specified gene and perform a model optimization.

orig_bounds = eo.gene_knock_outs('b0025')
solution = eo.optimize()
eo.set_variable_bounds(orig_bounds)
print(f'mutant gr: {solution.objective_value:.3f} h-1')
Parameters:

genes (str or list(str)) – gene identifiers

Returns:

original reaction bounds (lb, ub)

Return type:

dict(str, tuple)

moma(wt_fluxes, linear=False)[source]

Perform a MOMA based optimization using the gurobipy interface.

Supported for FBA and ECM models.

Ref: Segre et al. 2002, Analysis of optimality in natural and perturbed metabolic networks.

MOMA (linear = False):

minimize: ∑ (v - wt)^2

s.t.: S * v = 0; v_lb ≤ v ≤ v_ub

v: disturbed flux distribution; wt: wild type flux distribution.

This requires adding transformed variables vt with adjusted bounds, adding new constraints vt = v - wt and minimizing the objective function ∑ vt^2.

Linear MOMA (linear = True):

minimize: ∑ abs(v - wt)

s.t.: S * v = 0; v_lb ≤ v ≤ v_ub

introducing non-negative transformed flux variables and related constraints and minimizing ∑ (vtp + vtn)

vtp: [0 … (v_ub - wt)] with constraints vtp ≥ v - wt

vtn: [0 … (wt - v_lb)] with constraints vtn ≥ wt - v

Example: Perform MOMA analysis for a disturbed network, caused by gene deletion.

fo.medium = medium
wt_solution = fo.pfba()

orig_bounds = fo.gene_knock_outs('b0003')
moma_solution = fo.moma(wt_solution.fluxes)
fo.set_variable_bounds(orig_bounds)
print(f'{moma_solution.fluxes['BIOMASS_Ec_iML1515_core_75p37M']:.3f} h-1')
Parameters:

  • wt_fluxes (dict(str, float)) – wild type flux distribution, e.g. pFBA fluxes

  • linear (bool) – using quadratic (Euclidean distance) or linear formulation (Manhattan) (default: False)

Returns:

MOMA determined solution

Return type:

Solution

room(wt_fluxes, linear=False, delta=0.03, epsilon=0.001, time_limit=30.0)[source]

Perform a ROOM based optimization using the gurobipy interface.

Ref: Shlomi et al., 2005, Regulatory on off minimization of metabolic flux changes after genetic perturbations

ROOM (linear = False):

minimize: ∑ y^i

s.t. S * v = 0; v_lb ≤ v ≤ v_ub

v - y * (v_ub - wt_u) <= wt_u; wt_u = wt + delta * abs(wt) + epsilon

v - y * (v_lb - wt_l) <= wt_l; wt_l = wt - delta * abs(wt) - epsilon

v: disturbed flux distribution; wt: wild type flux distribution.

y_i ∈ {0,1}

Linear ROOM (linear = True):

0 <= y_i <= 1

Example: Perform ROOM analysis for a disturbed network, caused by gene deletion.

fo.medium = medium
wt_solution = fo.pfba()

orig_bounds = fo.gene_knock_outs('b0003')
room_solution = fo.room(wt_solution.fluxes)
fo.set_variable_bounds(orig_bounds)
print(f'{room_solution.fluxes['BIOMASS_Ec_iML1515_core_75p37M']:.3f} h-1')
Parameters:

  • wt_fluxes (dict(str, float)) – wild type flux distribution, e.g. pFBA fluxes

  • linear (bool) – using MILP (False) or relaxed LP (True) formulation (default: False)

  • delta (float) – relative tolerance range (default: 0.03)

  • epsilon (float) – absolute tolerance range (default: 1e-3)

  • time_limit (float) – time limit in seconds (default: 30.0)

Returns:

ROOM determined solution

Return type:

Solution

Solution

Solution object, aligned to COBRApy Solution.

class f2xba.optimization.optimize.Solution(status, objective_value=numpy.nan, fluxes=None, reduced_costs=None, shadow_prices=None)[source]

Structure holding information from optimization result.

Aligned to COBRApy cobra.Solution.

Example: Load and optimize a GECKO model using the gurobipy interface.

from f2xba import EcmOptimization

eo = EcmOptimization('iML1515_GECKO.xml')
eo.medium = {rid: 1000.0 for rid in lb_medium}
solution = eo.optimize()
solution.objective_value

Code for __repr__(), _repr_html_() copied from cobra.Solution.

status

Status string from optimization.

objective_value

Objective value (float).

fluxes

Values of optimization variables (pandas Series).

reduced_costs

Reduced cost of optimization variables (pandas Series).

shadow_prices

Shadow prices of constraints (pandas Series).