Reactions flipped and dropped from enzyme constraint

[VonAlphaBisZulu]

Hello,

Nice package, thank you very much!

There might be an issue with how you split your reactions that causes reversible parts of reactions to drop out from the enzyme constraint.

I was looking at the bounds of the enzyme-constraint reactions and noticed that cobrapy and optlang bounds are occasionally different. It seems like the reversible reactions get flipped once more when they go into optlang and then are no longer part of the enzyme constraint, so the model can use them in the reverse direction "for free". It's possible that I missed something but you may want to check.

import sys
sys.path.append('../../../../ECMpy/code/')
from cobrapy_ec_model_function import *
json_model_path="../../../../ECMpy/model/iML1515_irr_enz_constraint_adj_round2.json"
enz_model=get_enzyme_constraint_model(json_model_path)

# Look for reactions that are part of the enzyme constraint, but have contradictory optlang and cobrapy bounds
enzyme_constraint = [c for c in enz_model.constraints if c.ub>0][0]
reacs = []
for v in enzyme_constraint.variables:
    cp_ub = float(enz_model.reactions.get_by_id(v.name).upper_bound)
    cp_lb = float(enz_model.reactions.get_by_id(v.name).lower_bound)
    if v.ub != cp_ub or v.lb != cp_lb:
        reacs.append(v.name)
        print("Variable: "+v.name)
        print("optlang bounds: "+str(v))
        print("cobrapy bounds: "+str(cp_lb)+" <= "+str(v.name)+" <= "+str(cp_ub))
        if len(reacs) >= 3:
            break
        
print("=== variable lookup in enzyme constraint ===")

# Look if these reactions exist somewhere in the model under a different name
for r in reacs:
    for v in enz_model.variables:
        if r in str(v):
            print(("" if v in enzyme_constraint.variables else "NOT ")+"in enzyme constraint: "+str(v))

produces

Variable: SHK3Dr_reverse
optlang bounds: 0 <= SHK3Dr_reverse <= 0.0
cobrapy bounds: -1000.0 <= SHK3Dr_reverse <= 0.0
Variable: ACOAD2f
optlang bounds: 0 <= ACOAD2f <= 0.0
cobrapy bounds: -1000.0 <= ACOAD2f <= 0.0
Variable: ACACT4r_num2
optlang bounds: 0 <= ACACT4r_num2 <= 0.0
cobrapy bounds: -1000.0 <= ACACT4r_num2 <= 0.0
=== variable lookup in enzyme constraint ===
NOT in enzyme constraint: 0 <= SHK3Dr_reverse_d5c8f <= -0.0
in enzyme constraint: 0 <= SHK3Dr_reverse <= 0.0
NOT in enzyme constraint: 0 <= SHK3Dr_reverse_reverse_f4102 <= 1000.0
in enzyme constraint: 0 <= ACOAD2f <= 0.0
NOT in enzyme constraint: 0 <= ACOAD2f_reverse_6e942 <= 1000.0
in enzyme constraint: 0 <= ACACT4r_num2 <= 0.0
NOT in enzyme constraint: 0 <= ACACT4r_num2_reverse_47058 <= 1000.0

[VonAlphaBisZulu]

It might be easier to avoid the issues with cobrapy and optlang altogether if you express the enzyme constraint through pseudo-reactions and pseudo-variables. That way, everything is clean in the cobrapy and optlang will figure out the necessary details in the back.
The steps after splittling the reactions would be:

• Add a new metabolite (enzyme pool)
• Drain the new enzyme-pool metabolite in all enzyme-driven reactions, use your coefficients for that stoichiometry
• Add a pseudo reaction to produce/control the amount of available enzyme and bound that reaction with your enzyme-constraint-bound

[maozhitao]

Q1. There are occasional issues related to the disparities between cobrapy and optlang boundaries. When converting a cobrapy-formatted model to an ECMpy-formatted enzyme-constrained model, the first step involves splitting reversible reactions into two irreversible reactions. In this process, reactions with a flux range of (-1000, 1000) are considered reversible, while reactions with a flux range of (0, 1000) are treated as irreversible. However, due to the presence of reactions with a flux range of (-1000, 0) in the model, we mistakenly overlooked them, resulting in the aforementioned problem. Therefore, we have made corrections to the convert_to_irreversible(model) function, and after the code modification, this issue should no longer occur.

Q2: 2. Regarding the issue of random characters appearing in reaction names in variables: During the cobrapy calculation process, reversible reactions are split into two irreversible reactions for computation. In this splitting process, random characters (e.g., reverse_f4102) are automatically added after the reverse reaction to prevent naming conflicts. Since ECMpy models undergo the splitting of reversible reactions, the flux range for reverse reactions is (0, 0), and they are excluded from enzyme constraints. For metabolic network models, the same issue of adding random characters to reaction names can arise. Additionally, we experimented with the enzyme constraint formulation you suggested (pseudo-reactions and pseudo-variables). However, enzyme-constrained models generated using this approach also exhibit the problem of random characters in reaction names during computation. Therefore, this issue is inherent to the cobrapy model calculation mechanism rather than a problem with how you have modeled it.

[VonAlphaBisZulu]

The reverse reactions with random character sequences do not pose a problem. Since you split reactions into forward and reverse parts, most of the time, these reverse parts have boundaries (0,0) and are simply irrelevant. The main issue was (or is, I haven't checked, again), that there reactions whose reverse part should be limited by enzyme constraint but isn't.

The main reason I advocate for an implementation with pseudo metabolites and reactions is robustness and the improved universality and compatibility with plenty of other tools e.g., when you export sbml, when you use visualization like Escher or CNApy or compute intervention strategies with StrainDesign. The constraint-implementation that is currently used, is unique to Optlang and limited to vanilla cobrapy functions. If Optlang changes their syntax your package will break. If cobrapy uses another backend than Optlang, your package will break. You can see that it's even sensitive to how you split your reactions into forward and reverse. But metabolites, reactions and reaction boundaries are universal and will always be featured and supported. So my suggestion may help getting broadening your user base. E.g. someone might build an enzyme constraint model, export it as SBML and then load it in MATLAB for further computations.

[VonAlphaBisZulu]

Optlang constraints are conceptually and mathematically nice, but actually impractical

[VonAlphaBisZulu]

I checked your changes and they solved the issue! Well done and thanks :)

In case you're considering the pseudo reaction thing (I'd be happy if you did), take a look at pull request #5