2. Model preprocessing

A model must be thermodynamically sound in order to be used in PTA, which is usually not the case for published genome-scale models. We provide a simple method that automatically preprocesses a model for use in PTA.

pta.prepare_for_pta(model)

The modifications performed by this method include relaxing flux bounds, removing blocked reactions and removing “dead-end” futile cycles, thus they do not add further constraints to the capabilities of the model.

However, in some cases it may be relevant to obtain a complete report of the issues present in the model or to resolve these issues manually. The following sections explain how to do that. Only either prepare_for_pta() or the steps below are required.

2.1. Blocked reactions

In PTA all reactions have a well-defined direction, meaning that all reaction must have non-zero flux. This requires that all blocked reactions are removed from the model in advance. You can easily do this using COBRApy’s functions:

import cobra
from cobra.flux_analysis.variability import find_blocked_reactions

model.remove_reactions(
    find_blocked_reactions(model),
    remove_orphans=True
)

2.2. Forced internal cycles

Many models contain thermodynamically unfeasible internal cycles that must be active in any non-zero solution. This is usually a consequence of incorrect irreversibilities and “dead-end” cycles. For example:

The SUCDi and FRD7 reactions in the e_coli_core model both have the same stoichiometry and are irreversible, but in opposing directions.
```
SUCDi: succ_c + q8_c  -> q8h2_c + fum_c
FRD7:  q8h2_c + fum_c -> succ_c + q8_c
```
This would imply that both directions are thermodynamically favorable, which is not possible. The solution is to make both reactions reversible and let flux and thermodynamic constraints determine their direction.
The FOMETRi and THFAT in the iML1515 model suffer of the same problem:
```
THFAT:   h2o_c + methf_c -> 5fthf_c + h_c
FOMETRi: 5fthf_c + h_c   -> h2o_c + methf_c
```
However, this is a “dead-end” cycle, since only 5fthf_c does not participate to any other reaction. In this case, even making the two reactions reversible would not help, because any non-zero flux solution at steady state requires the two reactions to have opposing directions. Since cycles like this would anyway have no effect on the rest of the network we simply remove the reactions involved.

We emphasize that the presence of internal cycles in general is not a problem, and one of the strengths of PTA is to obtain loopless, thermodynamically feasible solutions. This section deals with internal cycles that must always be active and thus make the model unfeasible in principle.

2.2.1. Detecting inconsistencies

You can use the StructuralAssessment class to list all the problematic cycles in the model.

assessment = pta.StructuralAssessment(model, 'Biomass', 'ATPM')
assessment.summary()

> The following internal cycles are thermodynamically unfeasible, but must be active
in any non-zero flux solution, meaning that the model is thermodynamically
inconsistent.
...
0. THMDt2pp_copy1, THMDt2pp_copy2
1. ALAt2pp_copy1, ALAt2pp_copy2
2. CYTDt2pp_copy1, CYTDt2pp_copy2
3. ACCOAL, PPAKr, PTA2
4. ADNt2pp_copy1, ADNt2pp_copy2
5. FOMETRi, THFAT
...

Note that you need to specify the identifier of the biomass and ATP maintenance reactions.

2.2.2. Resolving inconsistencies

You can then resolve the problems manually or attempt to solve them automatically (logging is only needed to print the actions performed):

import logging
logging.basicConfig()
logging.getLogger().setLevel(logging.INFO)

assessment.autoresolve(model)

Removing bounds for all reactions in the cycle: THMDt2pp_copy1, THMDt2pp_copy2.
Removing bounds for all reactions in the cycle: ALAt2pp_copy1, ALAt2pp_copy2.
Removing bounds for all reactions in the cycle: CYTDt2pp_copy1, CYTDt2pp_copy2.
Removing bounds for all reactions in the cycle: ACCOAL, PPAKr, PTA2.
Removing bounds for all reactions in the cycle: ADNt2pp_copy1, ADNt2pp_copy2.
Removing all reactions of the dead-end cycle: FOMETRi, THFAT.
...

If you identified inconsistencies in a condition-specific model, you can apply autoresolve() on the base model in order to have the curation available for any condition-specific model generated from it.

2.3. Computing tighter bounds

To make the model simpler for later steps (PMO and TFS) it is good practice to restrict the bounds of each reaction to the bounds computed with Flux Variability Analysis.

pta.utils.tighten_model_bounds(model)

In order to not overconstrain the model and to avoid introducing small coefficients that would be challenging for PMO, this function may add small numerical margins. In any case, the resulting bounds are guaranteed to be equal or tighter than the original ones.