Merge pull request #34 from LCSB-BioCore/mk-convert-models

Properly support multiple model types and conversions

.gitignore

@@ -3,4 +3,6 @@
 docs/build
 test/data/ecoli_core.xml
 /vscode
-/.vscode
\ No newline at end of file
+/.vscode
+*.jl.*.cov
+.*.swp

src/base/solver.jl

@@ -30,7 +30,11 @@ end
 Use JuMP to optimize an instance of a [`MetabolicModel`](@ref). Returns a tuple
 that contains the new model and a vector of its variables.
 """
-function optimizeModel(model::LM, optimizer; sense = MOI.MIN_SENSE) where {LM<:MetabolicModel}
+function optimizeModel(
+    model::LM,
+    optimizer;
+    sense = MOI.MIN_SENSE,
+) where {LM<:MetabolicModel}
     optimization_model, x = makeOptimizationModel(model, optimizer; sense = sense)
     JuMP.optimize!(optimization_model)
     return (optimization_model, x)

src/base/utilities.jl

-function Base.isequal(model1::LinearModel, model2::LinearModel)
-    return model1.S == model2.S &&
-           model1.b == model2.b &&
-           model1.C == model2.C &&
-           model1.cl == model2.cl &&
-           model1.cu == model2.cu &&
-           model1.c == model2.c &&
-           model1.xl == model2.xl &&
-           model1.xu == model2.xu &&
-           model1.rxns == model2.rxns &&
-           model1.mets == model2.mets
-end
+Base.isequal(model1::LinearModel, model2::LinearModel) =
+    isequal(model1.S, model2.S) &&
+    isequal(model1.b, model2.b) &&
+    isequal(model1.c, model2.c) &&
+    isequal(model1.xl, model2.xl) &&
+    isequal(model1.xu, model2.xu) &&
+    isequal(model1.rxns, model2.rxns) &&
+    isequal(model1.mets, model2.mets)
 
-function Base.copy(model::LinearModel)
-    return LinearModel(
-        model.S,
-        model.b,
-        model.C,
-        model.cl,
-        model.cu,
-        model.c,
-        model.xl,
-        model.xu,
-        model.rxns,
-        model.mets,
-    )
-end
+Base.copy(model::LinearModel) =
+    LinearModel(model.S, model.b, model.c, model.xl, model.xu, model.rxns, model.mets)
+
+Base.isequal(model1::CoupledLinearModel, model2::CoupledLinearModel) =
+    isequal(model1.lm, model2.lm) &&
+    isequal(model1.C, model2.C) &&
+    isequal(model1.cl, model2.cl) &&
+    isequal(model1.cu, model2.cu)
+
+Base.copy(model::CoupledLinearModel) =
+    CoupledLinearModel(model.lm, model.C, model.cl, model.cu)

src/io/sbml.jl

@@ -2,35 +2,4 @@
 """
 Import a LinearModel from a SBML file
 """
-function loadSBMLModel(filename::String)::LinearModel
-    model = SBML.readSBML(filename)
-
-    mets, rxns, S = SBML.getS(model)
-    b = zeros(length(mets))
-    c = SBML.getOCs(model)
-
-    lbu = SBML.getLBs(model)
-    ubu = SBML.getUBs(model)
-
-    # SBML has units that we can't easily interpret yet, BUT the chances are
-    # that in a reasonable SBML file all units will be the same and we don't
-    # need to care about them at all. So we just check it here and throw an
-    # error if not.
-
-    unit = lbu[1][2]
-    getvalue = (val, _)::Tuple -> val
-    getunit = (_, unit)::Tuple -> unit
-
-    allunits = unique([getunit.(lbu) getunit.(ubu)])
-
-    if length(allunits) != 1
-        throw(
-            DomainError(
-                allunits,
-                "The SBML file uses multiple units; loading needs conversion",
-            ),
-        )
-    end
-
-    return LinearModel(S, b, c, getvalue.(lbu), getvalue.(ubu), rxns, mets)
-end
+loadSBMLModel(filename::String)::SBMLModel = SBMLModel(SBML.readSBML(filename))

src/reconstruction/coupling.jl

 """
-Verifies that vectors and matrices have the expected dimensions.
-"""
-function checkCouplingConstraintsInputDimensions(
-    m::LinearModel,
-    C::M,
-    cl::V,
-    cu::V,
-) where {M<:MatType,V<:VecType}
-
-    length(cu) == length(cl) ||
-        throw(DimensionMismatch("`cl` and `cu` don't have the same size"))
-    size(C) == (length(cl), nReactions(m)) || throw(
-        DimensionMismatch("`C` size should be length(cl) x the number of reactions in m"),
-    )
-end
-
-"""
-Add constraints of the following form to a LinearModel and return a modified one.
+Add constraints of the following form to a CoupledLinearModel and return a modified one.
 
 The arguments are same as for in-place `addCouplingConstraints!`.
 """
-function addCouplingConstraints(m::LinearModel, args...)
+function addCouplingConstraints(m::CoupledLinearModel, args...)
     newLp = deepcopy(m)
     addCouplingConstraints!(newLp, args...)
     return newLp
 end
 
+"""
+Add constraints to a plain `LinearModel` (converts it to the coupled type)
+"""
+addCouplingConstraints(m::LinearModel, args...) =
+    addCouplingConstraints(convert(CoupledLinearModel, m), args...)
+
 """
 In-place add coupling constraints in form
 ```
@@ -33,7 +22,7 @@ In-place add coupling constraints in form
 ```
 """
 function addCouplingConstraints!(
-    m::LinearModel,
+    m::CoupledLinearModel,
     c::V,
     cl::AbstractFloat,
     cu::AbstractFloat,
@@ -48,21 +37,20 @@ end
 
 
 function addCouplingConstraints!(
-    m::LinearModel,
+    m::CoupledLinearModel,
     C::M,
     cl::V,
     cu::V,
 ) where {M<:MatType,V<:VecType}
 
-    C = sparse(C)
-    cl = sparse(cl)
-    cu = sparse(cu)
-
-    checkCouplingConstraintsInputDimensions(m, C, cl, cu)
+    all([length(cu), length(cl)] .== size(C, 1)) ||
+        throw(DimensionMismatch("mismatched numbers of constraints"))
+    size(C, 2) == nReactions(m) ||
+        throw(DimensionMismatch("mismatched number of reactions"))
 
-    m.C = vcat(m.C, C)
-    m.cl = vcat(m.cl, cl)
-    m.cu = vcat(m.cu, cu)
+    m.C = vcat(m.C, sparse(C))
+    m.cl = vcat(m.cl, sparse(cl))
+    m.cu = vcat(m.cu, sparse(cu))
 end
 
 
@@ -71,7 +59,7 @@ Remove coupling constraints from the linear model and return the modified model.
 
 Arguments are the same as for in-place version `removeCouplingConstraints!`.
 """
-function removeCouplingConstraints(m::LinearModel, args...)
+function removeCouplingConstraints(m::CoupledLinearModel, args...)
     newModel = deepcopy(m)
     removeCouplingConstraints!(newModel, args...)
     return newModel
@@ -79,14 +67,14 @@ end
 
 
 """
-Removes a set of coupling constraints from a LinearModel in-place.
+Removes a set of coupling constraints from a CoupledLinearModel in-place.
 """
-function removeCouplingConstraints!(m::LinearModel, constraint::Int)
+function removeCouplingConstraints!(m::CoupledLinearModel, constraint::Int)
     removeCouplingConstraints!(m, [constraint])
 end
 
 
-function removeCouplingConstraints!(m::LinearModel, constraints::Vector{Int})
+function removeCouplingConstraints!(m::CoupledLinearModel, constraints::Vector{Int})
     toBeKept = filter(e -> e ∉ constraints, 1:nCouplingConstraints(m))
     m.C = m.C[toBeKept, :]
     m.cl = m.cl[toBeKept]
@@ -98,7 +86,7 @@ end
 Change the lower and/or upper bounds ('cl' and 'cu') for given coupling constraints
 """
 function changeCouplingBounds!(
-    model::LinearModel,
+    model::CoupledLinearModel,
     constraints::Vector{Int};
     cl::V = Array{Float64}(undef, 0),
     cu::V = Array{Float64}(undef, 0),

src/reconstruction/modeling.jl

@@ -104,12 +104,14 @@ function addReactions(
     xl = sparse(xl)
     xu = sparse(xu)
 
-    checkInputDimensions(Sp, b, c, xl, xu, rxns, mets)
+    all([length(b), length(mets)] .== size(Sp, 1)) ||
+        throw(DimensionMismatch("inconsistent number of metabolites"))
+    all(length.([c, xl, xu, rxns]) .== size(Sp, 2)) ||
+        throw(DimensionMismatch("inconsistent number of reactions"))
 
     newReactions = findall(Bool[!(rxn in m.rxns) for rxn in rxns])
     newMetabolites = findall(Bool[!(met in m.mets) for met in mets])
 
-
     if checkConsistency
         (newReactions1, newMetabolites1) =
             verifyConsistency(m, Sp, b, c, xl, xu, rxns, mets, newReactions, newMetabolites)
@@ -130,12 +132,13 @@ function addReactions(
 
     newS = hcat(extS, extSp)
     newb = vcat(m.b, b[newMetabolites])
-    newC = hcat(m.C, spzeros(size(m.C, 1), length(newReactions)))
+    #newC = hcat(m.C, spzeros(size(m.C, 1), length(newReactions)))
     newc = vcat(m.c, c[newReactions])
     newxl = vcat(m.xl, xl[newReactions])
     newxu = vcat(m.xu, xu[newReactions])
     newRxns = vcat(m.rxns, rxns[newReactions])
-    newLp = LinearModel(newS, newb, newC, m.cl, m.cu, newc, newxl, newxu, newRxns, newMets)
+    #newLp = LinearModel(newS, newb, newC, m.cl, m.cu, newc, newxl, newxu, newRxns, newMets)
+    newLp = LinearModel(newS, newb, newc, newxl, newxu, newRxns, newMets)
 
     if checkConsistency
         return (newLp, newReactions, newMetabolites)
@@ -144,56 +147,6 @@ function addReactions(
     end
 end
 
-"""
-Verifies that vectors and matrices have the expected dimensions.
-"""
-function checkInputDimensions(
-    S::M1,
-    b::V,
-    C::M2,
-    cl::V,
-    cu::V,
-    c::V,
-    xl::V,
-    xu::V,
-    rxns::K,
-    mets::K,
-) where {M1<:MatType,M2<:MatType,V<:VecType,K<:StringVecType}
-    n_c = length(c)
-
-    length(cu) == length(cl) ||
-        throw(DimensionMismatch("`cl` and `cu` don't have the same size"))
-    size(C) == (length(cu), n_c) ||
-        throw(DimensionMismatch("`C` shape doesn't match with `cu` and `c`"))
-
-    checkInputDimensions(S, b, c, xl, xu, rxns, mets)
-
-end
-
-function checkInputDimensions(
-    S::M,
-    b::V,
-    c::V,
-    xl::V,
-    xu::V,
-    rxns::K,
-    mets::K,
-) where {M<:MatType,V<:VecType,K<:StringVecType}
-
-    n_c = length(c)
-    n_b = length(b)
-
-    length(xu) == length(xl) ||
-        throw(DimensionMismatch("`xl` and `xu` don't have the same size"))
-    n_c == length(xl) || throw(DimensionMismatch("`c` doesn't have the same size as `xl`"))
-
-    size(S) == (n_b, n_c) ||
-        throw(DimensionMismatch("`S` shape doesn't match with `c` and `mets`"))
-
-    length(rxns) == n_c || throw(DimensionMismatch("`rxns` size doesn't match with `S`"))
-    length(mets) == n_b || throw(DimensionMismatch("`mets` size doesn't match with `S`"))
-end
-
 """
 Verifies the consistency of a given model
 """
@@ -250,14 +203,15 @@ function removeReactions(m::LinearModel, rxns::Vector{Int})
     unique!(metsToKeep)
     newS = m.S[metsToKeep, rxnsToKeep]
     newb = m.b[metsToKeep]
-    newC = m.C[:, rxnsToKeep]
+    #newC = m.C[:, rxnsToKeep]
     newc = m.c[rxnsToKeep]
     newxl = m.xl[rxnsToKeep]
     newxu = m.xu[rxnsToKeep]
     newRxns = m.rxns[rxnsToKeep]
     newMets = m.mets[metsToKeep]
     newModel =
-        LinearModel(newS, newb, newC, m.cl, m.cu, newc, newxl, newxu, newRxns, newMets)
+    #LinearModel(newS, newb, newC, m.cl, m.cu, newc, newxl, newxu, newRxns, newMets)
+        LinearModel(newS, newb, newc, newxl, newxu, newRxns, newMets)
     return newModel
 end
 
@@ -271,7 +225,7 @@ function removeReactions(m::LinearModel, rxn::String)
 end
 
 
-function removeReactions(m::LinearModel, rxns::Array{String,1})
+function removeReactions(m::LinearModel, rxns::Vector{String})
     rxnIndices = [findfirst(isequal(name), m.rxns) for name in intersect(rxns, m.rxns)]
     if isempty(rxnIndices)
         return m

src/types/linearModel.jl

@@ -6,16 +6,12 @@ A concrete linear optimization problem of the form:
 ```
 min c^T x
 s.t. S x = b
-    cₗ ≤ C x ≤ cᵤ
-    xₗ ≤ x ≤ xᵤ
+      xₗ ≤ x ≤ xᵤ
 ```
 """
 mutable struct LinearModel <: MetabolicModel
     S::SparseMat
     b::SparseVec
-    C::SparseMat
-    cl::SparseVec
-    cu::SparseVec
     c::SparseVec
     xl::SparseVec
     xu::SparseVec
@@ -32,43 +28,13 @@ mutable struct LinearModel <: MetabolicModel
         mets::K,
     ) where {V<:VecType,M<:MatType,K<:StringVecType}
 
-        sS = sparse(S)
-        sb = sparse(b)
-        sC = spzeros(0, length(c))
-        scl = spzeros(0)
-        scu = spzeros(0)
-        sc = sparse(c)
-        sxl = sparse(xl)
-        sxu = sparse(xu)
+        all([length(b), length(mets)] .== size(S, 1)) ||
+            throw(DimensionMismatch("inconsistent number of reactions"))
 
-        LinearModel(sS, sb, sC, scl, scu, sc, sxl, sxu, rxns, mets)
-    end
-
-    function LinearModel(
-        S::M1,
-        b::V,
-        C::M2,
-        cl::V,
-        cu::V,
-        c::V,
-        xl::V,
-        xu::V,
-        rxns::K,
-        mets::K,
-    ) where {V<:VecType,M1<:MatType,M2<:MatType,K<:StringVecType}
+        all([length(c), length(xl), length(xu), length(rxns)] .== size(S, 2)) ||
+            throw(DimensionMismatch("inconsistent number of reactions"))
 
-        checkInputDimensions(S, b, C, cl, cu, c, xl, xu, rxns, mets)
-
-        sS = sparse(S)
-        sb = sparse(b)
-        sC = sparse(C)
-        scl = sparse(cl)
-        scu = sparse(cu)
-        sc = sparse(c)
-        sxl = sparse(xl)
-        sxu = sparse(xu)
-
-        new(sS, sb, sC, scl, scu, sc, sxl, sxu, rxns, mets)
+        new(sparse(S), sparse(b), sparse(c), sparse(xl), sparse(xu), rxns, mets)
     end
 end
 
@@ -129,17 +95,36 @@ end
 """
     coupling(a::LinearModel)::SparseMat
 
-Coupling constraint matrix for a `LinearModel`.
+Coupling constraint matrix for a `LinearModel`, actually empty.
 """
 function coupling(a::LinearModel)::SparseMat
-    a.C
+    spzeros(0, nReactions(a))
 end
 
 """
     couplingBounds(a::LinearModel)::Tuple{SparseVec,SparseVec}
 
-Coupling bounds for a `LinearModel`.
+Coupling bounds for a `LinearModel`, in fact always empty.
 """
 function couplingBounds(a::LinearModel)::Tuple{SparseVec,SparseVec}
-    (a.cl, a.cu)
+    (spzeros(0), spzeros(0))
+end
+
+
+"""
+    Base.convert(::Type{LinearModel}, m::M) where {M <: MetabolicModel}
+
+Make a `LinearModel` out of any compatible model type.
+"""
+function Base.convert(::Type{LinearModel}, m::M) where {M<:MetabolicModel}
+    (xl, xu) = bounds(m)
+    LinearModel(
+        stoichiometry(m),
+        balance(m),
+        objective(m),
+        xl,
+        xu,
+        reactions(m),
+        metabolites(m),
+    )
 end

src/types/linearModelCoupled.jl

+
+"""
+    struct CoupledLinearModel <: MetabolicModel
+
+The linear model with additional coupling constraints in the form
+```
+    cₗ ≤ C x ≤ cᵤ
+```
+"""
+mutable struct CoupledLinearModel <: MetabolicModel
+    lm::LinearModel
+    C::SparseMat
+    cl::SparseVec
+    cu::SparseVec
+
+    function CoupledLinearModel(
+        lm::MetabolicModel,
+        C::M,
+        cl::V,
+        cu::V,
+    ) where {V<:VecType,M<:MatType}
+
+        length(cu) == length(cl) ||
+            throw(DimensionMismatch("`cl` and `cu` need to have the same size"))
+        size(C) == (length(cu), nReactions(lm)) ||
+            throw(DimensionMismatch("wrong dimensions of `C`"))
+
+        new(convert(LinearModel, lm), sparse(C), sparse(cl), sparse(cu))
+    end
+end
+
+reactions(a::CoupledLinearModel) = reactions(a.lm)
+metabolites(a::CoupledLinearModel) = metabolites(a.lm)
+stoichiometry(a::CoupledLinearModel) = stoichiometry(a.lm)
+bounds(a::CoupledLinearModel) = bounds(a.lm)
+balance(a::CoupledLinearModel) = balance(a.lm)
+objective(a::CoupledLinearModel) = objective(a.lm)
+
+"""
+    coupling(a::CoupledLinearModel)::SparseMat
+
+Coupling constraint matrix for a `CoupledLinearModel`.
+"""
+function coupling(a::CoupledLinearModel)::SparseMat
+    a.C
+end
+
+"""
+    nCouplingConstraints(a::CoupledLinearModel)::Int
+
+The number of coupling constraints in a `CoupledLinearModel`.
+"""
+function nCouplingConstraints(a::CoupledLinearModel)::Int
+    return size(a.C, 1)
+end
+
+"""
+    couplingBounds(a::CoupledLinearModel)::Tuple{SparseVec,SparseVec}
+
+Coupling bounds for a `CoupledLinearModel`.
+"""
+function couplingBounds(a::CoupledLinearModel)::Tuple{SparseVec,SparseVec}
+    (a.cl, a.cu)
+end
+
+"""
+    Base.convert(::Type{CoupledLinearModel}, m::M) where {M <: MetabolicModel}
+
+Make a `CoupledLinearModel` out of any compatible model type.
+"""
+function Base.convert(::Type{CoupledLinearModel}, m::M) where {M<:MetabolicModel}
+    (cl, cu) = couplingBounds(m)
+    CoupledLinearModel(convert(LinearModel, m), coupling(m), cl, cu)
+end

src/types/sbmlModel.jl

+"""
+    struct SBMLModel
+
+Thin wrapper around the model from SBML.jl library. Allows easy conversion from
+SBML to any other model format.
+"""
+struct SBMLModel <: MetabolicModel
+    m::SBML.Model
+end
+
+reactions(a::SBMLModel)::Vector{String} = [k for k in keys(a.m.reactions)]
+metabolites(a::SBMLModel)::Vector{String} = [k for k in keys(a.m.species)]
+nReactions(a::SBMLModel)::Int = length(a.m.reactions)
+nMetabolites(a::SBMLModel)::Int = length(a.m.species)
+
+"""
+    stoichiometry(a::SBMLModel)::SparseMat
+
+Recreate the stoichiometry matrix from the SBML model.
+"""
+function stoichiometry(a::SBMLModel)::SparseMat
+    _, _, S = SBML.getS(a.m)
+    return S
+end
+
+"""
+    bounds(a::SBMLModel)::Tuple{SparseVec,SparseVec}
+
+Get the lower and upper flux bounds of a `SBMLModel`. Throws `DomainError` in
+case if the SBML contains mismatching units.
+"""
+function bounds(a::SBMLModel)::Tuple{SparseVec,SparseVec}
+    lbu = SBML.getLBs(a.m)
+    ubu = SBML.getUBs(a.m)
+
+    unit = lbu[1][2]
+    getvalue = (val, _)::Tuple -> val
+    getunit = (_, unit)::Tuple -> unit
+
+    allunits = unique([getunit.(lbu) getunit.(ubu)])
+    length(allunits) == 1 || throw(
+        DomainError(
+            allunits,
+            "The SBML file uses multiple units; loading needs conversion",
+        ),
+    )
+
+    return sparse.((getvalue.(lbu), getvalue.(ubu)))
+end
+
+balance(a::SBMLModel)::SparseVec = spzeros(nMetabolites(a))
+objective(a::SBMLModel)::SparseVec = SBML.getOCs(a.m)
+coupling(a::SBMLModel)::SparseMat = spzeros(0, nReactions(a))
+nCouplingConstraints(a::SBMLModel)::Int = 0
+couplingBounds(a::SBMLModel)::Tuple{SparseVec,SparseVec} = (spzeros(0), spzeros(0))

test/base/utilities.jl

@@ -12,4 +12,5 @@
 
     cp = test_coupledLP()
     @test nCouplingConstraints(cp) == 2000
+    @test isequal(cp, copy(cp))
 end

test/data/testModels.jl

@@ -41,15 +41,17 @@ test_sparseLP() = LinearModel(
     ["m$x" for x = 1:4000],
 )
 
-test_coupledLP() = LinearModel(
-    sprand(4000, 3000, 0.5),
-    sprand(4000, 0.5),
+test_coupledLP() = CoupledLinearModel(
+    LinearModel(
+        sprand(4000, 3000, 0.5),
+        sprand(4000, 0.5),
+        sprand(3000, 0.5),
+        sprand(3000, 0.5),
+        sprand(3000, 0.5),
+        ["r$x" for x = 1:3000],
+        ["m$x" for x = 1:4000],
+    ),
     sprand(2000, 3000, 0.5),
     sprand(2000, 0.5),
     sprand(2000, 0.5),
-    sprand(3000, 0.5),
-    sprand(3000, 0.5),
-    sprand(3000, 0.5),
-    ["r$x" for x = 1:3000],
-    ["m$x" for x = 1:4000],
 )

test/io/sbml.jl

@@ -12,15 +12,20 @@ if cksum != "78692f8509fb36534f4f9b6ade23b23552044f3ecd8b48d84d484636922ae907"
     @warn "The downloaded E Coli core model seems to be different from the expected one. Tests may fail." cksum
 end
 
-@testset "SBML import" begin
-    m = loadSBMLModel(sbmlfile)
+@testset "SBML import and conversion" begin
+    sbmlm = loadSBMLModel(sbmlfile)
+    m = convert(LinearModel, sbmlm)
 
-    @test size(m.C) == (0, 95)
-    @test size(m.S) == (92, 95)
+    @test size(stoichiometry(sbmlm)) == (92, 95)
+    @test size(stoichiometry(m)) == (nMetabolites(sbmlm), nReactions(sbmlm))
     @test length(m.S.nzval) == 380
-    @test length(m.b) == 92
+    @test length.(bounds(sbmlm)) == (95, 95)
+    @test length.(bounds(m)) == (95, 95)
     @test all([length(m.xl), length(m.xu), length(m.c)] .== 95)
 
-    @test m.mets[1:3] == ["M_succoa_c", "M_ac_c", "M_fru_b"]
-    @test m.rxns[1:3] == ["R_EX_fum_e", "R_ACONTb", "R_GLNS"]
+    @test metabolites(m)[1:3] == ["M_succoa_c", "M_ac_c", "M_fru_b"]
+    @test reactions(m)[1:3] == ["R_EX_fum_e", "R_ACONTb", "R_GLNS"]
+
+    cm = convert(CoupledLinearModel, sbmlm)
+    @test nCouplingConstraints(cm) == 0
 end

test/reconstruction/coupling.jl

 @testset "Coupling constraints" begin
-    cp = test_LP()
-    @test size(cp.S) == (4, 3)
-    newCp = addCouplingConstraints(cp, cp.S[end, :], -1.0, 1.0)
+    cp = convert(CoupledLinearModel, test_LP())
+    @test size(cp.lm.S) == (4, 3)
+    @test size(stoichiometry(convert(LinearModel, cp))) == (4, 3)
+    newCp = addCouplingConstraints(cp, stoichiometry(cp)[end, :], -1.0, 1.0)
     @test nCouplingConstraints(cp) + 1 == nCouplingConstraints(newCp)
 
-    newCp = addCouplingConstraints(cp, cp.S[1:2, :], [-1.0; -1.0], [1.0; 1.0])
+    newCp = addCouplingConstraints(cp, stoichiometry(cp)[1:2, :], [-1.0; -1.0], [1.0; 1.0])
     @test nCouplingConstraints(cp) + 2 == nCouplingConstraints(newCp)
 
     nC = nCouplingConstraints(cp)
-    addCouplingConstraints!(cp, cp.S[end, :], -1.0, 1.0)
+    addCouplingConstraints!(cp, stoichiometry(cp)[end, :], -1.0, 1.0)
     @test nC + 1 == nCouplingConstraints(cp)
-    addCouplingConstraints!(cp, cp.S[1:2, :], [-1.0; -1.0], [1.0; 1.0])
+    addCouplingConstraints!(cp, stoichiometry(cp)[1:2, :], [-1.0; -1.0], [1.0; 1.0])
     @test nC + 3 == nCouplingConstraints(cp)
 
     nC = nCouplingConstraints(cp)
@@ -23,6 +24,7 @@
     cp = test_coupledLP()
     nC = nCouplingConstraints(cp)
     newCp = removeCouplingConstraints(cp, 1)
+    @test size(coupling(cp)) == (nC, nReactions(cp))
     @test nC - 1 == nCouplingConstraints(newCp)
     @test nCouplingConstraints(cp) == nC
     newCp = removeCouplingConstraints(cp, [1, 2])
@@ -33,8 +35,10 @@
 
     cp = test_coupledLP()
     changeCouplingBounds!(cp, [3, 1], cl = [-10.0, -20], cu = [10.0, 20])
-    @test cp.cl[[1, 3]] == [-20, -10]
-    @test cp.cu[[1, 3]] == [20, 10]
+    cl, cu = couplingBounds(cp)
+    @test cl[[1, 3]] == [-20, -10]
+    @test cu[[1, 3]] == [20, 10]
     changeCouplingBounds!(cp, [1000, 1001], cl = [-50.0, -60.0])
-    @test cp.cl[[1000, 1001]] == [-50.0, -60.0]
+    cl, cu = couplingBounds(cp)
+    @test cl[[1000, 1001]] == [-50.0, -60.0]
 end