from IPython.display import display
import re

from cameo import models
from cameo.strain_design import pathway_prediction

model = models.bigg.iMM904

predictor = pathway_prediction.PathwayPredictor(model=model, compartment_regexp=re.compile(".*_c$"))

pathways = predictor.run(product="vanillin", max_predictions=5)
Pathway 1
id equation lower_bound upper_bound
0 MNXR230 H2O + 3,4-dihydroxybenzoate + NADP(+) <=> O2 +... -1000 1000
1 MNXR5336 vanillin + H2O + NAD(+) <=> 2.0 H(+) + vanilla... -1000 1000
2 MNXR5340 H2O + NAD(+) + 3,4-dihydroxybenzoate + formald... -1000 1000
Pathway 2
id equation lower_bound upper_bound
0 MNXR5336 vanillin + H2O + NAD(+) <=> 2.0 H(+) + vanilla... -1000 1000
1 MNXR5340 H2O + NAD(+) + 3,4-dihydroxybenzoate + formald... -1000 1000
2 MNXR68718 H2O + 3,4-dihydroxybenzoate <=> 3-dehydroshiki... -1000 1000
Pathway 3
id equation lower_bound upper_bound
0 MNXR5336 vanillin + H2O + NAD(+) <=> 2.0 H(+) + vanilla... -1000 1000
1 MNXR5340 H2O + NAD(+) + 3,4-dihydroxybenzoate + formald... -1000 1000
2 MNXR14769 H2O + 3,4-dihydroxybenzoate + NAD(+) <=> O2 + ... -1000 1000
Pathway 4
id equation lower_bound upper_bound
0 MNXR5336 vanillin + H2O + NAD(+) <=> 2.0 H(+) + vanilla... -1000 1000
1 MNXR5340 H2O + NAD(+) + 3,4-dihydroxybenzoate + formald... -1000 1000
2 MNXR5836 3.0 H(+) + anthranilate + O2 + NADPH <=> CO(2)... -1000 1000
3 MNXR7067 H(+) + 3,4-dihydroxybenzoate <=> CO(2) + catechol -1000 1000
Pathway 5
id equation lower_bound upper_bound
0 MNXR5336 vanillin + H2O + NAD(+) <=> 2.0 H(+) + vanilla... -1000 1000
1 MNXR5340 H2O + NAD(+) + 3,4-dihydroxybenzoate + formald... -1000 1000
2 MNXR6101 3.0 H(+) + anthranilate + O2 + NADH(2-) <=> CO... -1000 1000
3 MNXR7067 H(+) + 3,4-dihydroxybenzoate <=> CO(2) + catechol -1000 1000 
from cameo import phenotypic_phase_plane
from cameo.util import TimeMachine
from cameo.visualization.plotting import Grid

with Grid(nrows=3) as grid:
    for i, pathway in enumerate(pathways):
        with TimeMachine() as tm:
            pathway.plug_model(model, tm=tm)
            ppp = phenotypic_phase_plane(model, variables=[model.reactions.biomass_SC5_notrace], objective=pathway.product)
            ppp.plot(grid=grid, width=450, height=350, title="Pathway %i" % (i+1), axis_font_size="12pt")