Source-linked AI summary
Creation and analysis of biochemical constraint-based models: the COBRA Toolbox v3.0
Laurent Heirendt, Sylvain Arreckx, Thomas Pfau, Sebastián N. Mendoza, Anne Richelle, Almut Heinken, Hulda S. Haraldsdóttir, Jacek Wachowiak, Sarah M. Keating, Vanja Vlasov, Stefania Magnusdóttir, Chiam Yu Ng, German Preciat, Alise Žagare, Siu H. J. Chan, Maike K. Aurich, Catherine M. Clancy, Jennifer Modamio, John T. Sauls, Alberto Noronha, Aarash Bordbar, Benjamin Cousins, Diana C. El Assal, Luis V. Valcarcel, Iñigo Apaolaza, Susan Ghaderi, Masoud Ahookhosh, Marouen Ben Guebila, Andrejs Kostromins, Nicolas Sompairac, Hoai M. Le, Ding Ma, Yuekai Sun, Lin Wang, James T. Yurkovich, Miguel A. P. Oliveira, Phan T. Vuong, Lemmer P. El Assal, Inna Kuperstein, Andrei Zinovyev, H. Scott Hinton, William A. Bryant, Francisco J. Aragón Artacho, Francisco J. Planes, Egils Stalidzans, Alejandro Maass, Santosh Vempala, Michael Hucka, Michael A. Saunders, Costas D. Maranas, Nathan E. Lewis, Thomas Sauter, Bernhard Ø. Palsson, Ines Thiele, Ronan M. T. Fleming
TL;DR
Genome-scale molecular measurements have expanded faster than the ability to derive mechanistic insight from them, while linear approximations are insufficient for some biochemical processes. The COBRA Toolbox 3.0 updates and expands interoperable constraint-based modelling methods, data integration, and numerical optimisation capabilities. It supports model generation, analysis, and data integration across diverse molecular systems biology applications, while requiring users to evaluate candidate reactions and reproduce analyses through code.
Problem
Genome-scale molecular measurements generate abundant data, but deriving comprehensive mechanistic insight remains challenging, and linear approximations are insufficient for some biochemical processes.
Method
The protocol updates COBRA Toolbox versions 1.0 and 2.0 with expanded constraint-based modelling methods, optimisation solvers, model formats, computational functions, and metabolomic-data integration.
Results
The protocol produces COBRA model structures and files and provides candidate reactions for completing metabolic reconstructions.
Takeaways & Limitations
COBRA Toolbox 3.0 offers a flexible suite of constraint-based modelling methods that users can combine for tailored analyses across biochemical networks.
Takeaways & Limitations
The COBRA Toolbox focuses on constraint-based reconstruction and analysis rather than general-purpose computational biology, and graphically driven tools do not provide computational specifications for reproducing point-and-click analyses.
Abstract
from arXiv · showhide
COnstraint-Based Reconstruction and Analysis (COBRA) provides a molecular mechanistic framework for integrative analysis of experimental data and quantitative prediction of physicochemically and biochemically feasible phenotypic states. The COBRA Toolbox is a comprehensive software suite of interoperable COBRA methods. It has found widespread applications in biology, biomedicine, and biotechnology because its functions can be flexibly combined to implement tailored COBRA protocols for any biochemical network. Version 3.0 includes new methods for quality controlled reconstruction, modelling, topological analysis, strain and experimental design, network visualisation as well as network integration of chemoinformatic, metabolomic, transcriptomic, proteomic, and thermochemical data. New multi-lingual code integration also enables an expansion in COBRA application scope via high-precision, high-performance, and nonlinear numerical optimisation solvers for multi-scale, multi-cellular and reaction kinetic modelling, respectively. This protocol can be adapted for the generation and analysis of a constraint-based model in a wide variety of molecular systems biology scenarios. This protocol is an update to the COBRA Toolbox 1.0 and 2.0. The COBRA Toolbox 3.0 provides an unparalleled depth of constraint-based reconstruction and analysis methods.
INTRODUCTION
The COBRA Toolbox 3.0 addresses the challenge of extracting mechanistic insight from expanding molecular datasets by providing interoperable constraint-based reconstruction and analysis methods. It updates earlier Toolbox versions with broader modelling functionality, narrative tutorials, open-source development, and support for diverse biochemical modelling applications.
- Genome-scale molecular measurements have expanded rapidly, but deriving biological insight from integrated experimental data remains challenging.
- Mechanistic COBRA models incorporate prior biochemical-network information to connect genotype and phenotype, including when mechanistic knowledge is incomplete.
- COBRA Toolbox 3.0 updates versions 1.0 and 2.0 with new standards, solvers, modelling methods, data integration, thermodynamic constraints, visualisation, and kinetic modelling.
- Narrative tutorials, open-source collaboration, and community support facilitate reuse and adaptation of COBRA methods into customised protocols.
- The Toolbox provides interoperable methods across reconstruction, analysis, design, visualisation, data integration, and supporting numerical functions through one software environment.
- The Toolbox is focused on constraint-based reconstruction and analysis rather than general-purpose computational biology, and some published methods remain unavailable as direct MATLAB implementations.
MATERIALS
The COBRA Toolbox requires a MATLAB-compatible computing environment, the Toolbox itself, shell and version-control utilities, and at least one compatible optimization solver. Additional software and MATLAB toolboxes are needed only for specific optional functions, while solver and MATLAB-version compatibility can constrain use.
- A 64-bit Intel or AMD computer needs at least 8 GB of RAM and 10 GB of free hard-drive storage.
- Installation and contribution require a stable internet connection, the COBRA Toolbox 3.0 or newer, a bash or UNIX shell, curl 7.0 or newer, and git 1.8 or newer.Windows requires Git Bash integration, while macOS requires Xcode.
- The required platform includes Linux, macOS, or Windows compatible with MATLAB R2014b or newer, with R2016a or newer needed for live scripts.Tutorials can also run on R2014b using provided .m files.
- MATLAB releases and solver binaries may temporarily lack compatibility, and some solvers have additional system requirements or platform restrictions.The dqqMinos and Minos solvers may only be used on Unix, and some COBRA.binary executables may not work on a given system.
- Optional analyses may require libSBML, MATLAB toolboxes, Chemaxon Calculator Plugins, Java, Python, NumPy, OpenBabel, or Reaction Decoder Tool.These dependencies support model exchange, topology, flux variability analysis, sampling, cheminformatics, and related functionality.
- At least one linear-programming solver is required for basic constraint-based modelling, while GLPK and several other solvers are installed by default.On Windows, the OPTI solver suite must be installed separately to use its interface.
PROCEDURE
The procedure initializes the COBRA Toolbox, verifies dependencies and solvers, tests functionality, and imports, validates, exports, and extends biochemical models. It supports multiple model formats and rBioNet-based reaction addition, while consistency checks remain essential for legacy or malformed files.
- Initialization: Each MATLAB session begins by initializing the COBRA Toolbox, either automatically for regular use or manually when contributing or managing multiple repositories.Initialization checks software dependencies and reports them in the command window.
- Initialization: Initialization selects default optimization solvers, using Gurobi when installed and otherwise GLPK for LP and MILP problems.Users should inspect the reported solver compatibility and availability.
- Testing: The optional test suite verifies core-function execution and outputs on the locally configured system.
- Model input and validation: Models can be read from SBML, Excel, Simpheny, and MATLAB files, with readCbModel checking consistency and converting legacy model structures.Using readCbModel instead of load supports future-proof model handling.
- Model input and validation: Legacy MATLAB models may contain deprecated or invalid fields that readCbModel cannot always correct, requiring verifyModel and manual adjustment.After correction, the model should be saved again and loaded with readCbModel.
- Model input and validation: After reading a model, the MATLAB workspace should contain a struct with core fields including S, lb, ub, c, osense, b, csense, rxns, mets, genes, and rules.
- Model output: Models can be exported with writeCbModel to SBML or MATLAB files, provided the structure and field data satisfy the required specifications.The output file is placed in the indicated folder or current directory; verifyModel is recommended when writing fails.
- Reconstruction editing: rBioNet supports adding or removing reactions and gene-reaction associations through its graphical interface or database commands, followed by conversion to the current model structure.The relevant metabolite, reaction, and compartment databases must be available and correctly configured.
22 | Flux balance analysis,
This section presents COBRA procedures for computing flux-balance solutions, handling infeasible constraints, identifying sparse pathways and network gaps, and selecting candidate gap-filling reactions. It emphasizes that optimal flux vectors may be non-unique and that proposed gap fills require biological evaluation.
- optimizeCbModel computes flux balance analysis and related variants through changes to optional input arguments.
- The optimal objective value is unique under feasible constraints, but the optimal flux vector is usually not unique and should not be interpreted in isolation.A biologically interpreted flux vector should be unique to an optimization problem.
- A strictly convex regularization term can produce a unique optimal flux vector while preserving the original objective.Setting minNorm to 10^-6 implements a quadratic regularization, which may need adjustment in high-dimensional models.
- Relaxed flux balance analysis addresses infeasible constraints by minimizing bound relaxations, with priorities set according to biochemical context.Relaxations should be minimized and justified for each adjusted bound.
- Sparse flux balance analysis minimizes the number of active reactions subject to stoichiometric, bound, and optionally optimal-objective constraints.The resulting vector represents a minimal-support steady-state pathway under the specified conditions.
- Gap-filling methods identify dead-end metabolites and blocked reactions, producing candidate reactions that must be assessed for biological and physiological plausibility.Dead-end metabolites mark network gaps, while candidate additions require evaluation before inclusion.
Extracellular metabolomic data
The COBRA Toolbox integrates extracellular and intracellular metabolomic measurements into condition-specific or dynamic metabolic models. These workflows constrain models using media composition, measured exchange fluxes, time-course concentration changes, and selected omics-supported reactions.
- Extracellular metabolomic data: MetaboTools integrates extracellular metabolomic data with metabolic models to generate experimentally testable mechanistic hypotheses.The integration seeks to include measured metabolites while adding minimal uptake and secretion metabolites.
- Extracellular metabolomic data: Medium composition is encoded as external reaction bounds, including water and bicarbonate needed to support biomass synthesis.The resulting modelMedium represents the defined fresh medium.
- Extracellular metabolomic data: Metabolomic matrices represent uptake as negative flux and secretion as positive flux, with consistent flux units across the model.Rows correspond to exchanged metabolites and columns correspond to samples.
- Extracellular metabolomic data: setQuantConstraints generates functional, condition-specific models from medium constraints, metabolomic measurements, tolerances, growth requirements, and secretion or uptake restrictions.Each sample receives a model whose constraints reflect its medium and measured extracellular metabolomic data.
- Intracellular metabolomic data: uFBA integrates absolutely quantified time-course metabolomics into metabolic models, including intracellular data collected when metabolite levels are not at steady state.Its workflow estimates concentration-change slopes and integrates them as intracellular flux constraints.
- Intracellular metabolomic data: buildUFBAmodel uses measured metabolites, concentration-change slopes, confidence intervals, and ignored slopes to produce a constrained COBRA model and identify deviations from steady state.The output records integrated metabolites and relaxed nodes, and may add sink or demand reactions to retain feasibility.
- Integration of transcriptomic and proteomic data: Context-specific reconstruction methods extract subsets of generic models using experimental data, with algorithms differing in how they select or retain reactions.FASTCORE supports core reactions, GIMME minimizes low-expression usage, and related methods implement alternative expression-based trade-offs.
Adding biological constraints to a flux balance model
The protocol adds experimentally and biologically motivated constraints to metabolic models, then tests whether constrained models reproduce expected biochemical functions and quantitative behavior.
- Cell-type models can become condition-specific by imposing biomass, environmental, or maximal enzyme-activity constraints.
- Literature-derived bounds must use the model’s flux units, and curating biomass composition and precursor turnover rates can take 4–8 weeks.Once constraints are available, imposing the corresponding reaction bounds takes less than 5 minutes.
- A degradation constraint can require the summed flux through selected reactions to meet a specified lower bound.The example represents phosphatidylcholine degradation through PCHOLP_hs, PLA2_2, and SMS.
- Fidelity testing combines stoichiometric checks, biochemical tests, and comparisons between model predictions and independent experimental data.The protocol tests whether closed models cannot produce ATP from water and whether expected metabolic functions remain feasible.
- Quantitative validation compares predicted ATP yields from different carbon sources and oxygen conditions with biochemical-literature values.
MinSpan Pathways: a sparse basis of the nullspace of a stoichiometric matrix
The COBRA Toolbox represents feasible biochemical behavior through pathway and variability analyses, including sparse nullspace bases, flux ranges, and uniform sampling of high-dimensional solution spaces.
- MinSpan pathways: MinSpan computes a sparse set of n−rank(S) linearly independent flux modes forming a basis of the stoichiometric right nullspace.It uses a greedy sequence of mixed-integer linear optimization problems.
- Flux variability analysis: Flux balance analysis generally admits infinitely many steady-state flux vectors with the same optimal objective, motivating flux variability analysis.
- Flux variability analysis: The Toolbox provides standard and advanced flux variability interfaces that return minimum and maximum reaction fluxes under configurable optimality conditions.The standard output is a pair of n-dimensional minimum and maximum flux vectors.
- High-dimensional analysis: High-dimensional variability analysis can use fastFVA with CPLEX or distributedFBA.jl across parallel processors or cluster workers.The distributed workflow supports solvers including GLPK, CPLEX, Clp, Gurobi, and Mosek.
- Sampling: Uniform sampling characterizes feasible steady-state flux vectors, with parameter tuning recommended for higher-dimensional models.Models with fewer than 2000 variables can usually use default sampling settings.
- Sampling: CHRR rounds anisotropic solution spaces before coordinate hit-and-run sampling, producing more distributed and smoother converged reaction distributions.The example uses nSkip = 8 × (dim(fluxSpace))^2 and nSamples = 1000.
- Strain design: The Toolbox implements strain-design algorithms including OptKnock, OptGene, GDLS, and OptForce.OptForce additionally identifies up- and down-regulation strategies alongside gene deletions.
Atomically resolve a metabolic reconstruction
Atomic resolution links metabolic stoichiometry to molecular structures, enabling atom mapping, conserved-moiety analysis, and thermodynamic constraints on reaction feasibility.
- Genome-scale stoichiometric matrices encode molecular species and reactions without explicitly representing metabolite atomic structures.
- Metabolite structures can be obtained from databases, cheminformatics tools, or manual drawing before atom mapping.
- Atom mappings can be predicted computationally, retrieved from curated databases, or drawn manually from reaction stoichiometry and metabolite structures.RDT is described as suitable for genome-scale mapping, while CLCA handles explicit hydrogen atoms in fully protonated reactions.
- Canonicalized atom mappings support conserved-moiety identification, where invariant molecular substructures are tracked across network transformations.The workflow constructs an atom transition network and identifies moieties from its structure.
- Thermochemical integration estimates transformed Gibbs energies using metabolite structures, training data, compartment conditions, and concentration bounds.Multi-compartment models must be processed as entire networks so intercompartment thermodynamic potentials are included.
- Thermodynamic constraints assign reaction directionality and can disallow stoichiometrically balanced loops, yielding thermodynamically feasible steady-state fluxes.The resulting solutions satisfy energy conservation and the second law of thermodynamics.
Convert a flux balance model into a kinetic model
The protocol motivates kinetic modeling because linear flux balance analysis can be insufficient for nonlinear biochemical processes, then states structural conditions for internally consistent kinetic models.
- Flux balance analysis uses linear optimization to analyze genome-scale biochemical networks and can produce biologically meaningful predictions in some situations.
- Variational kinetic modeling is motivated for biochemical processes where a linear approximation is insufficient.
- The kinetic formulation assumes a stoichiometrically consistent network with forward and reverse stoichiometric matrices representing reversible reactions.The stated conditions include m < n and sparse stoichiometric matrices.
Compute a non-equilibrium kinetic steady state
The protocol computes non-equilibrium kinetic steady states by solving nonlinear systems derived from elementary forward and reverse reaction kinetics. The COBRA Toolbox uses solver-based workflows that accept a kinetic model, solver, initial point, and algorithm parameters, while supporting steady-state visualisation through ReconMap.
- Kinetic steady-state formulation: Elementary forward and reverse reaction rates are represented as exponential functions of kinetic parameters and molecular-species concentrations.The formulation uses component-wise exp and ln functions applied to kinetic parameters and concentration vectors.
- Kinetic steady-state formulation: A concentration vector c∗ is a steady state if and only if it satisfies f(c∗) = 0.
- Numerical solution methods: Sparsity in stoichiometric matrices and non-unique local zeros motivate specialised algorithms beyond generic nonlinear least-squares methods.The section discusses derivative-free methods and adaptive Levenberg–Marquardt approaches, with line-search or trust-region globalisation for arbitrary starting points.
- Numerical solution methods: The optimizeVKmodels workflow requires a vKModel containing F and R, a nonlinear-system solver, an initial point x0, and solver parameters.Optional settings include the maximum number of iterations, adaptive mode, and kinetic parameters; defaults are used otherwise.
- Visualisation: ReconMap visualisation places context-specific flux predictions in biochemical pathway maps, helping interpret model predictions within their biological context.The workflow uses MINERVA to create overlays, which can then be viewed in a browser and exported as PDF or PNG graphics.
Contributing to the COBRA Toolbox with MATLAB.devTools
MATLAB.devTools provides a branch-and-fork workflow for contributing code to the COBRA Toolbox while reducing the risk that new changes break heavily used functions. Contributors work from the develop branch, publish changes to personal forks, and use pull requests for acceptance into openCOBRA.
- Workflow rationale: Continuous integration linked to git is used because frequent, interlinked contributions can break heavily used code.
- Repository setup: Users work from the stable master branch, whereas contributors create branches from develop in a personal fork.The openCOBRA repository is read-only for contributors, while the fork provides read and write access.
- Repository setup: The contribute procedure clones or updates the fork, checks out develop, and initializes a new contribution branch for code changes.The fork is created in a local fork-cobratoolbox folder, and contributors can then add, modify, or delete files.
- Updating contributions: Procedure retrieves changes from openCOBRA and rebases an existing contribution after the develop branch.Contributors should publish their commits before continuing work on an existing feature.
- Publishing contributions: Publishing uploads commits to the personal fork; acceptance into the official repository requires a pull request.Contributors can publish without immediately opening a pull request, then open one later when the contribution is ready.
Engaging with the COBRA Toolbox forum
The COBRA Toolbox community supports users through documentation, an online forum, and direct contributions. The forum is especially useful for finding solutions to recurring problems and recent installation or configuration issues.
- Finding help: The FAQ section is a starting point for answers to questions and issues encountered with the COBRA Toolbox.
- Finding help: The public forum enables users to search for solutions to previously recognised problems similar to their own.It is particularly useful for recent installation and configuration issues caused by asynchronous development of integrated software packages.
- Posting questions: Forum questions should be detailed and include the output of generateSystemConfigReport to help maintainers understand the system configuration.Membership application is required before posting to reduce spam.
- Community participation: Community contributions are welcomed, and responses can generally be expected within 1-2 days when posting guidelines are followed.
? TROUBLESHOOTING
The troubleshooting guidance maps common COBRA Toolbox problems to likely causes and corrective actions. It emphasises checking installation and configuration, validating model data, adjusting numerical settings, and using community support when needed.
- Installation and initialization: Initialization warnings may result from binaries incompatible with the operating system, so software requirements and supported versions should be checked first.The documentation recommends reading initialization output literally and testing the toolbox after verifying installations.
- Model and file validation: Model validation should precede export because missing required fields or invalid data can affect reconstructions and models.The verifyModel(model) function reports invalid data, while accepted file formats and annotations have specific compatibility conditions.
- Model consistency: Incorrect reaction identification or balance checks can arise from unexpected reaction prefixes or missing and incorrectly specified formulae.The troubleshooting steps recommend inspecting reaction formulae and rerunning consistency analysis after omitting obviously imbalanced reactions.
- Numerical troubleshooting: Sampling problems may reflect short runtimes or low nSkip and nSamples settings, which should be increased until marginal flux distributions are smooth and unimodal.
- Development tools: Contribution procedures can fail when forks are unreachable, SSH keys are misconfigured, local forks are outdated, or local changes conflict with repository updates.Recommended remedies include checking connectivity and keys, updating or recloning forks, preserving backups, and resetting after failed rebases when appropriate.
TIMING
The protocol reports execution times for COBRA Toolbox workflows spanning reconstruction, model testing, optimisation, design, thermodynamic and kinetic modelling, visualisation, and development.
- Reconstruction: Adding reactions with rBioNet takes 1–10^3 s, whereas spreadsheet-based addition takes 1–10^3s.
- Model quality control: Scaling checks, precision-solver selection, consistency analyses, and detection of leaks, dead ends, and blocked reactions take 1–10^5 s.
- Model testing: Biological constraint addition and qualitative or quantitative fidelity testing require approximately 10^2–10^3 s.
- Flux analysis: Flux analyses range from 1–10^3 s for low-dimensional variability analysis to 1–10^5 s for high-dimensional analysis.
- Design and mechanistic modelling: Strain-design searches and atomically resolved reconstruction require 10–10^5 s, while thermodynamic and kinetic modelling steps generally require 1–10^3 s.
- Visualisation and development: Network visualisation takes 1–10^3 s for Paint4Net and 1–10^2 s for ReconMap, while MATLAB.devTools contributions take 1–30 s.
ANTICIPATED RESULTS
The anticipated results illustrate Toolbox setup, solver configuration, model loading and checking, biochemical fidelity testing, flux analysis, strain design, and network visualisation. Examples include ATP-yield agreement, identification of an energy-generating cycle, and OptForce interventions.
- Initialisation: Initialisation checks dependencies, updates submodules, configures compatible solvers, and reports solver status and update options.
- Solver configuration: The Toolbox reports solver support across LP, MILP, QP, MIQP, and NLP classes, with Gurobi supporting all listed classes except NLP.
- Model handling: The workflow loads reconstructions or models into a MATLAB structure named model and can export the model information to a specified file.
- Numerical precision: Double-precision and quad-precision workflows are selected through scaling estimates, while dqqMinos provides feasibility and optimality tolerance of 10^-15.
- Analytical examples: The examples identify an energy-generating cycle, an intuitive succinate transporter solution, and OptForce interventions for aerobic succinate overproduction in E. coli.
- Biochemical fidelity testing: 32 ATP per unit of glucose is predicted for Recon3Dmodel, compared with 31 ATP from biochemical literature.
Supplementary information
The supplementary information provides introductory guidance for MATLAB, terminal or shell use, and contributing to the COBRA Toolbox.
- The supplementary material includes a brief tutorial on MATLAB software.
- It also provides guidelines for using the terminal or shell.
- Additional guidance covers using git and contributing to the COBRA Toolbox.
Supplementary Manual 1 - MATLAB basics
The supplementary MATLAB manual introduces basic command-line, file-management, matrix, vector, plotting, and statistical operations used in Toolbox workflows.
- MATLAB basics: The manual introduces MATLAB startup, function syntax, command-line notation, and script execution.
- File and workspace management: It demonstrates loading and saving workspaces and variables, clearing the workspace, and deleting individual variables.
- Data structures and linear algebra: Matrix and vector examples cover constructing matrices, creating word-entry vectors, calculating matrix rank, and computing a null space.
- Visualisation: The manual also introduces singular value decomposition, graphics windows, plotting, histograms, axis labels, titles, and multiple graphs.
Supplementary Manual 2 - Shell or Terminal basics
The manual introduces shell commands for navigating, managing, and synchronizing a forked COBRA Toolbox repository. It then outlines branch-based contribution, rebasing, cherry-picking, and file-history workflows.
- Shell or Terminal basics: Shell commands support directory navigation, file management, and retrieval of the current working path.The manual covers ls, cd, pwd, mkdir, rm -rf, and rm.
- Repository setup: Users clone their fork, initialize submodules, configure the upstream remote, and fetch upstream changes.These steps prepare the local fork and connect it to the main repository.
- Repository synchronization: The master and develop branches are updated by merging upstream branches and pushing the results to the fork.The manual gives separate command sequences for both branches.
- Contributions: Contributions are made on a new branch, where changes are reviewed, staged selectively, committed, pushed, and submitted as a pull request.The manual warns against using git add . because it can include hidden files and binaries.
- Branch maintenance: Rebasing incorporates upstream changes into an existing branch, while conflicts may require a merge tool and force-pushing afterward.The workflow requires checking for uncommitted changes before rebasing.
- Selective changes and history: Cherry-pick transfers a selected commit between branches, and file history can be inspected through git or MATLAB.devTools.The manual uses SHA1-based cherry-picking and provides git log --follow and history commands.