How do I apply all experimental flux measurements without over-constraining my model?

[zealottt]

Hey guys,

I am trying to set up a constraint-based model, question is the number of fluxes measured is large which means it is an over-determined system, I want to know how to cope with situation like this without expanding my model size?

[Midnighter]

Just a comment: You will need to provide a lot more details for us to be able to say something helpful/meaningful. In most cases, genome-scale models are under-determined that is why an optimization of an objective function approach is used. You are saying that you have more measured fluxes than reactions in your system? What is the model for in that case? How many measured fluxes and reactions do you have?

[zealottt]

Sorry, I should say my problem in details. I have 41 metabolites and 51 reactions in my model, but I can test 14 reaction flux, which makes it an over-determined system. Is there a way to take use of these 14 reactions as constraints without dumping any one of them?

[Midnighter]

1. The first thing that I would check is the relationship between uptake rate(s) and growth rate, for example, in a minimal medium with glucose as the sole carbon source, does the model predict a growth rate similar to the measured one or are there large differences? If there are large differences, you may have to adjust either the uptake or the growth rate as seems sensible (best discuss this with the experimentalist who measured the fluxes).
2. Typically, flux measurements are performed in triplicate. You could try setting the bounds of the reactions to mean flux - two standard deviations and the upper bound to mean flux + two standard deviations. If that still over-constrains your system, the last thing I would try is:
3. Constrain the uptake and growth rates. Then formulate an optimization problem that (given those constraints) minimizes the distance between a feasible flux distribution and the measured fluxes. I have some code for this so let me know if you need it.

[zealottt]

I am interested in the third point you put forward. What I can't understand is how to get a feasible flux in an over-determined system, something like randomly choose 10 of 14 reactions to get the flux vector and then pick up one among the vectors which has the least distance from the measured flux? I am curious about the code if you can share me with. Thank you.

[Midnighter]

I'm somewhat hesitant to make this code public because it is not comprehensively tested and also does not work under all circumstances. Please evaluate the results that you get with this carefully. I highlight this part of the license

WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND

adjust_fluxes2model is the function that you're looking for.

# -*- coding: utf-8 -*-

# Copyright 2017 Novo Nordisk Foundation Center for Biosustainability,
# Technical University of Denmark.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""Minimize the distance between experimental and simulated fluxes."""

from __future__ import absolute_import, division

import logging

from numpy import isinf
from pandas import Series
from cobra.util import add_absolute_expression
from optlang.symbolics import Zero

__all__ = ("constrain_with_observations", "minimize_distance",
           "adjust_fluxes2model")

LOGGER = logging.getLogger(__name__)


def constrain_with_observations(model, observations, index=None, lower=None,
                                upper=None, mean=None, variance=None,
                                inplace=False):
    """
    Constrain a metabolic model with observed fluxes.
    Will either use a lower and upper value to set the flux bounds or a mean
    value plus/minus some arbitrary measure of variance. Either the ``lower``
    or ``mean`` are required arguments.
    Parameters
    ----------
    model : cobra.Model
        The constraint-based metabolic model.
    observations : pandas.DataFrame
        Needs to contain at least one column with observed fluxes.
    index : str, optional
        The column name containing reaction identifiers. By default uses the
        index of the data frame.
    lower : str, optional
        The column name of the lower observed flux bounds.
    upper : str, optional
        The column name of the upper observed flux bounds.
    mean : str, optional
        The column name of the mean observed fluxes.
    variance : str, optional
        The column name of the variance around the mean observed fluxes.
    inplace : bool, optional
        Whether to return a constrained copy (default) or adjust the given
        model.
    Returns
    -------
    cobra.Model or None
        A copy of the original model constrained with the observations.
    """
    index, lower, upper = lower_upper_bounds(observations, index, lower, upper,
                                             mean, variance)
    if inplace:
        constr_model = model
    else:
        constr_model = model.copy()
    for rxn_id in index:
        try:
            rxn = constr_model.reactions.get_by_id(rxn_id)
        except KeyError:
            LOGGER.warning("'%s' not found in the model. Ignored.", rxn_id)
            continue
        rxn.bounds = lower[rxn_id], upper[rxn_id]
    if not inplace:
        return constr_model


def minimize_distance(model, observations, linear=True):
    """
    Find a flux distribution within a minimal distance to observed fluxes.
    Uses either a quadratic formulation (Manhattan distance) where the sum of
    the squares of the differences between simulated flux values and the
    actually observed fluxes is minimized or a linearized formulation where the
    sum of the absolute values of the differences between simulated flux values
    and the actually observed fluxes is minimized.
    Parameters
    ----------
    model : cobra.Model
        The constraint-based metabolic model.
    observations : pandas.Series
        A measured flux distribution (or any kind of distribution) whose
        index should contain reaction identifiers.
    linear : bool, optional
        Whether or not to use a linear formulation. If ``False`` requires a
        solver capable of optimizing quadratic problems.
    Returns
    -------
    cobra.Solution
        The result of the distance minimization.
    """
    prob = model.problem
    v = prob.Variable("moma_old_objective")
    c = prob.Constraint(model.solver.objective.expression - v, lb=0.0, ub=0.0,
                        name="old_objective_constraint")
    to_add = [v, c]
    new_obj = Zero
    with model:
        for rxn_id, flux in observations.iteritems():
            try:
                rxn = model.reactions.get_by_id(rxn_id)
            except KeyError:
                LOGGER.warning("'%s' not found in the model. Ignored.", rxn_id)
                continue
            if linear:
                components = add_absolute_expression(
                    model, rxn.flux_expression, name="moma_dist_" + rxn_id,
                    difference=flux, add=False)
                to_add.extend(components)
                new_obj += components.variable
            else:
                dist = prob.Variable("moma_dist_" + rxn.id)
                const = prob.Constraint(rxn.flux_expression - dist,
                                        lb=flux, ub=flux,
                                        name="moma_constraint_" + rxn.id)
                to_add.extend([dist, const])
                new_obj += dist**2
        model.add_cons_vars(to_add)
        model.objective = prob.Objective(new_obj, direction="min")
        solution = model.optimize(raise_error=True)
    return solution


def adjust_fluxes2model(model, observations, uncertainties=None, linear=True,
                        big_m=1E05):
    """
    Minimize the distance to observed fluxes accounting for multiple directions.
    If your observations include uncertainties the objective function, i.e.,
    minimizing the distance to the observations, is weighted by the inverse
    of the uncertainties.
    Parameters
    ----------
    model : cobra.Model
        The metabolic model
    observations : pandas.Series
        The observed fluxes. The index should contain reaction identifiers.
    uncertainties : pandas.Series, optional
        The uncertainties of the individual measurements, e.g., standard
        error. The index of the series should correspond at least partially
        to the ``observations``.
    linear : bool, optional
        Whether to minimize the linear or quadratic distance.
    big_m : float, optional
        Big M method value. This is used to resolve greater than inequalities
        and should be an adequately large number.
    Returns
    -------
    cobra.Solution
    """
    flux_col = "flux"
    weight_col = "weight"
    if uncertainties is None:
        data = observations.to_frame().join(Series([], name=weight_col))
    else:
        uncertainties.name = weight_col
        data = observations.to_frame().join(uncertainties)
    data.columns = [flux_col, weight_col]
    # replace missing and zero values
    data.loc[data[weight_col].isnull() | isinf(data[weight_col]) |
             (data[weight_col] == 0), weight_col] = 1
    prob = model.problem
    to_add = list()
    new_obj = list()
    with model:
        model.objective = prob.Objective(prob.Zero, direction="min", sloppy=True)
        for rxn_id, flux, weight in data[[flux_col, weight_col]].itertuples():
            try:
                rxn = model.reactions.get_by_id(rxn_id)
            except KeyError:
                LOGGER.warning("'%s' not found in the model. Ignored.", rxn_id)
                continue
            direction = prob.Variable("direction_" + rxn_id, type="binary")
            dist = prob.Variable("dist_" + rxn_id)
            forward_pos = prob.Constraint(
                flux - rxn.flux_expression - big_m * (1 - direction) - dist,
                ub=0, name="forward_pos_" + rxn_id)
            forward_neg = prob.Constraint(
                rxn.flux_expression - flux - big_m * (1 - direction) - dist,
                ub=0, name="forward_neg_" + rxn_id)
            reverse_pos = prob.Constraint(
                (-flux) - rxn.flux_expression - big_m * direction - dist,
                ub=0, name="reverse_pos_" + rxn_id)
            reverse_neg = prob.Constraint(
                rxn.flux_expression - (-flux) - big_m * direction - dist,
                ub=0, name="reverse_neg_" + rxn_id)
            if linear:
                new_obj.append((dist, 1 / weight))
            else:
                raise NotImplementedError("Need to set coefficients.")
#                new_obj.append((dist * dist, 1 / weight))
            to_add.extend([direction, dist, forward_pos, forward_neg,
                           reverse_pos, reverse_neg])
        model.add_cons_vars(to_add)
        model.objective.set_linear_coefficients({v: f for v, f in new_obj})
        solution = model.optimize(raise_error=True)
    return solution


def lower_upper_bounds(data, index=None, lower=None, upper=None, mean=None,
                       variance=None):
    if lower is not None:
        lower = data[lower]
        if upper is not None:
            upper = data[upper]
        else:
            upper = lower
    elif mean is not None:
        if variance is not None:
            lower = data[mean] - data[variance]
            upper = data[mean] + data[variance]
        else:
            lower = data[mean]
            upper = lower
    else:
        raise ValueError(
            "Either the lower flux value or the mean flux are required.")
    if index is None:
        index = data.index
    else:
        index = data[index]
    return index, lower, upper

[zealottt]

Thank you for your sharing！I will have a try but it may take some time.