This is package aims to provide constraint based reconstruction and analysis tools at the exa-scale in Julia.

## Installation

To install this package: `] add ???`. See the documentation for more information.

## Quick Example

Let's use `COBREXA.jl` to perform classic flux balance analysis on an *E. coli* community.

```julia

usingCOBREXA

usingJuMP

usingTulip# pick any solver supported by JuMP

# Import E. coli models (models have pretty printing)

model_1=read_model("iJO1366.json")

model_2=read_model("iJO1366.json")

model_3=read_model("iJO1366.json")

# Build an exascale model

exascale_model=join(model_1,model_2,model_3,...)

```

More funcionality is described in the documention, e.g. model construction and analysis in pure Julia.

### Citations

1) Ebrahim, A., Lerman, J.A., Palsson, B.O. & Hyduke, D. R. (2013). COBRApy: COnstraints-Based Reconstruction and Analysis for Python. BMC Systems Biology, 7(74). https://doi.org/10.1186/1752-0509-7-74

2) Heirendt, L., Arreckx, S., Pfau, T. et al. (2019). Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. Nat Protoc 14, 639–702. https://doi.org/10.1038/s41596-018-0098-2

3) Noor, E., Bar-Even, A., Flamholz, A., Lubling, Y., Davidi, D., & Milo, R. (2012). An integrated open framework for thermodynamics of reactions that combines accuracy and coverage. Bioinformatics, 28(15), 2037–2044. https://doi.org/10.1093/bioinformatics/bts317

4) Chang, A., Jeske, L., Ulbrich, S., Hofmann, J., Koblitz, J., Schomburg, I., Neumann-Schaal, M., Jahn, D., Schomburg, D.. (2021). BRENDA, the ELIXIR core data resource in 2021: new developments and updates. Nucleic Acids Research, 49(D1). https://doi.org/10.1093/nar/gkaa1025

A selection of standard COBRA functions have been implemented to make basic model analysis more convenient.

Additionally, `CobraTools.jl` allows you to easily formulate your own optimization problems using the structure of a constraint based model.

This makes it easy to experiment with custom algorithms etc.

## Flux balance analysis (FBA)

Flux balance analysis solves the linear program,

```math

\begin{aligned}

& \underset{v}{\text{max}}

& \sum_{i}^I {w_i}{v_i} \\

& \text{s. t.}

& Sv = 0 \\

& & v_{\text{LB}} \leq v \leq v_{\text{UB}} \\

\end{aligned}

```

using any `JuMP` compatible solver. Typically ``I`` is a singleton set that only includes the index of the biomass objective function, ``v_\mu``, and weight, ``w_\mu=1``.

```@docs

fba

```

Here, we use `Tulip.jl`, a pure Julia interior point linear program solver, with the `fba` function from `CobraTools.jl`.

Note that the result of the optimization problem, `sol` maps fluxes to reaction `id`s in the model, to simplify down stream analysis.

Since the fluxes are mapped to a dictionary, it makes it simple to export them as a JSON file for visualization in, e.g. [Escher](https://escher.github.io/#/).

```@example fba

using JSON

open("fluxes.json", "w") do io

JSON.print(io, sol)

end

```

## Solution inspection

Sometimes it is useful to investigate which reactions consume or produce a certain metabolite, or which exchange reactions are active.

This functionality is exposed via `metabolite_fluxes` and `exchange_reactions`.

Parsimonious FBA (pFBA) solves a two stage optimization problem. First, a classic FBA problem is solved to identify the unique maximum of the objective.

However, it should be noted that the fluxes from FBA are not unique (i.e. many fluxes may yield the objective optimum).

To yield a unique set of fluxes, and remove internal futile cycles, a secondary quadratic problem is imposed *ad hoc*.

Suppose that FBA has found the optimum of ``v_\mu = \mu``, pFBA then solves,

```math

\begin{aligned}

& \underset{v}{\text{min}}

& \sum_{i} {v_i^2} \\

& \text{s. t.}

& Sv = 0 \\

& & v_{\text{LB}} \leq v \leq v_{\text{UB}} \\

& & v_\mu = \mu

\end{aligned}

```

again using any JuMP compatible solver(s). In the `CobraTools.jl` implementation of pFBA, both the FBA and QP problem are solved internally in `pfba`, using similar input arguments as in `fba`. If multiple solvers are given, the first solver is used to solve the LP, and the second solver the QP, otherwise the same solver is used to solve both problems.

This is useful if the QP solver does not handle the LP problem well, as with OSQP.

An alternative, related formulation of this idea exists, called "CycleFreeFlux".

This replaces the quadratic formulation with an L1 (taxi cab) norm. While this should also remove futile cycles, it doesn't have the same uniqueness qualities and doesn't

really have much benefit beyond only solving two linear programs. See [Building your own optimization analysis script](@ref) for ideas about how to implement this yourself if you really need it.

```@docs

pfba

```

Here, we use `Tulip.jl` followed by `OSQP.jl`, with the `pfba` function from `CobraTools.jl`. Note that `OSQP.jl` has iffy performance, and is only included here because it is open source. We recommend that a commercial solver, e.g. `Gubobi.jl`, be used to simplify your user experience.

```@example fba

using OSQP

atts = Dict("eps_abs" => 5e-4,"eps_rel" => 5e-4, "max_iter" => 100_000, "verbose"=>false) # set solver attributes for QSQP

Flux variability analysis can also be used to investigate the degeneracy associated with flux balance analysis derived solutions (see also [Sampling Tools](@ref)).

`CobraTools.jl` exposes `fva` that sequentially maximizes and minimizes each reaction in a model subject to the constraint that each optimization problem also satisfies an initial FBA type objective optimum, below denoted by ``v_{\mu}=\mu``,

```math

\begin{aligned}

& \underset{v}{\text{max or min}}

& v_i \\

& \text{s. t.}

& Sv = 0 \\

& & v_{\text{LB}} \leq v \leq v_{\text{UB}} \\

& & v_{\mu} = \mu \\

\end{aligned}

```

```@docs

fva

```

## Building your own optimization analysis script

`CobraTools.jl` also makes it simple to construct customized optimization problems by making judicious use of [JuMP](https://jump.dev/).

Convenience functions make optimization problem construction, modification and data extraction from JuMP result objects easy.

Numerous external databases exist that can add functionality to constraint based models.

Currently, `CobraTools.jl` tries to make it easier to incorporate two such databases in analyses.

The first is Brenda, a database that maps enzymes to kinetic data.

The second is Equilibrator, a database that allows one to calculate thermodynamic information, e.g. ``\Delta G``, from reactions.

`CobraTools.jl` includes lightweight interfaces to these sources.

## Brenda Interface

To use the Brenda interface functions, you will need to download the database as a txt file [available here](https://www.brenda-enzymes.org/download_brenda_without_registration.php)(~250 MB).

Once the database has been downloaded, it can be parsed by `parse_brenda`.

```@docs

parse_brenda

```

This function returns an array of `BrendaEntry` structs, which are composed of `EnzymeParams` for each field extracted from the Brenda database.

Currently, only the ID (=EC number), TN (=turn over number), KM (=Michaelis-Menten constant, ``K_M``), KI (=Inhibition term for Michaelis-Menten kintics), and KKM (=ratio of TN/KM) numbers for each enzyme class (ID, or EC number) are extracted. All the structs have pretty printing enabled.

The Equilibrator interface requires that the Equilibrator-API has been installed and can be accessed through Julia's PyCall package. Refer to the [Equilibrator-API website](https://gitlab.com/equilibrator/equilibrator-api) for installation instructions. Within Julia, if you can call `pyimport("equilibrator_api")` successfully, then you will be able to use the functions exposed here. To actually use the functions insert `using PyCall` in your main level script (before or after `using CobraTools`).

@@ -22,4 +22,50 @@ If you want to contribute, please read these guidelines first:

```@contents

Pages = ["howToContribute.md"]

```

\ No newline at end of file

```

## Contents

```@contents

Pages = [

"model_structure.md",

"io.md",

"model_construction.md",

"basic_analysis.md",

"sampling_tools.md",

"external_tools.md",

"thermodynamics.md"

]

Depth=2

```

## Installation

To install this package: `] add https://github.com/stelmo/CobraTools.jl`.

Some of the optional features used in this package require external programs and/or data to be available. These are described below:

* The Equilibrator interface requires that the Equilibrator-API has been installed and can be accessed through Julia's PyCall package. Refer to the [Equilibrator-API website](https://gitlab.com/equilibrator/equilibrator-api) for installation instructions. Within Julia, if you can call `pyimport("equilibrator_api")` successfully, then you will be able to use the functions exposed here. To actually use the functions insert `using PyCall` in your main level script (before or after `using CobraTools`).

* To extract turnover numbers, Km, Kcat/Km and Ki from the Brenda database, you will need to download the database as a txt file [available here](https://www.brenda-enzymes.org/download_brenda_without_registration.php)(~250 MB).

The optimization solvers are implemented through `JuMP` and thus this package should be solver agnostic. All tests are conducted using `Tulip.jl` and `OSQP.jl`, but other solvers should also work (I mostly use `Gurobi.jl`).

## Quick Example

Let's perform flux balance analysis on a constraint based model.

adp = Metabolite("adp") # define another metabolite

metdict = Dict(atp => -1.0, adp => 1.0) # nb stoichiometries need to be floats

rxn = Reaction("dummy rxn", metdict, "for")

rxn.annotation["ec-code"] = ["0.0.0.0"]

rxn.grr = [[gene]] # only gene1 is required for this reaction to work

rxn # pretty printing

```

See the discussion in [Model Structure](@ref) about how to assign `grr` to a reaction.

Yet another way of defining a reaction is through overloading of the operators: `*, +, ∅, ⟶, →, ←, ⟵, ↔, ⟷`.

The longer and shorter arrows mean the same thing, i.e. `⟶` is the same as `→`, etc.

The other fields of the reaction still need to be set directly.

```@example

using CobraTools # hide

atp = Metabolite("atp") # hide

atp.name = "Adenosine triphosphate" # hide

atp.formula = "C10H12N5O13P3" # hide

atp.charge = -4 # hide

adp = Metabolite("adp") # hide

adp.formula = "C10H12N5O10P2" # hide

another_rxn = 2.0adp ⟶ 2.0*atp # forward reaction

another_rxn.id = "another dummy rxn"

another_rxn

```

When building exchange, demand, or sinks reactions the `∅` empty metabolite should be used to indicate that a metabolite is being created or destroyed.

```@example

using CobraTools # hide

adp = Metabolite("adp") # hide

adp.formula = "C10H12N5O10P2" # hide

ex_rxn = ∅ ⟷ adp # exchange reaction

```

It is also possible to check if a reaction is mass balanced by using `is_mass_balanced(rxn::Reaction)`.

Note, this function requires that all the metabolites in the reaction have formulas assigned to them to work properly.