#!/usr/bin/env python
# coding: utf-8

# # Optimize co-factor swap

# Many metabolic enzymes depend on co-factors to function. Keeping balance between co-factors is important for homeostasis and that balance might interact unfavorable with metabolic engineering. And example of such balance is that between the two similar co-factor pairs NAD+/NADH and NADP+/NADPH. These co-factors are not only similar there are even enzymes that catalyze the same reaction but depend on different co-factors. Enzymes can also be enigneered to change their co-factor preference.
# 
# This opens an opportunity for using co-factor swaps to optimize production of a target metabolite. Figuring out which reactions should be subjected to co-factor swap can be done using the [OptSwap algorithm](http://www.sciencedirect.com/science/article/pii/S1096717614000688). Briefly, the algorithm uses a genetic algorithm to test combinations of reactions to co-factor swap and then reports those that results in higher theoretical maximum yield.
# 
# We have implemented a variant of OptSwap in cameo and here we test it to reproduce the results in that paper for the aerobic yeast metabolism.
# 
# 

# Get the model and import modules.

# In[2]:


from cameo import models
model_orig = models.bigg.iJO1366


# In[3]:


from cameo.strain_design.heuristic.evolutionary.optimization import CofactorSwapOptimization
from cameo.strain_design.heuristic.evolutionary.objective_functions import product_yield
from cameo.strain_design.heuristic.evolutionary.objective_functions import biomass_product_coupled_yield
from cameo.util import TimeMachine
from cameo.flux_analysis.analysis import flux_variability_analysis as fva


# We make a copy of the model for easy testing.

# In[4]:


model = model_orig.copy()


# Make model changes as indicated in the paper.

# In[6]:


for rid in ['FHL', 'CAT', 'SPODM', 'SPODMpp']:
    model.reactions.get_by_id(rid).knock_out()
model.reactions.POR5.lower_bound = 0
model.reactions.EX_glc__D_e.lower_bound = -10
model.reactions.EX_o2_e.lower_bound = -10


# In[7]:


model.reactions.BIOMASS_Ec_iJO1366_core_53p95M.lower_bound = 0.1


# In the paper they get 0.77 maximum product yield for l-threonine, which we also get

# In[8]:


model.objective = model.reactions.EX_thr__L_e
(model.solve().f * 4) / (model.reactions.EX_glc__D_e.flux * 6)


# Let's run optswap using the `CofactorSwapOptimization` class in cameo. We use a product-yield function to evaluate how good a solution is, but there are other possibilities like `biomass_coupled_product_yield`.
# 
# GAPD is suggested in the paper as a suitable reactions

# In[9]:


py = product_yield(model.reactions.EX_thr__L_e, model.reactions.EX_glc__D_e)
optswap = CofactorSwapOptimization(model=model, objective_function=py)


# We also observe GAPD among the best options even though we considered many more reactions than in the original paper

# In[10]:


optswap.run(max_evaluations=2000, max_size=2)


# The created `optswap` class has properties to check which reactions were tested and to perform swapping, let's list the first 10.

# In[16]:


list(optswap.model.swapped_reactions)[0:10]


# In[32]:


optswap.model.reactions.EX_thr__L_e.model = optswap.model
optswap.model.objective = optswap.model.reactions.EX_thr__L_e
original = (optswap.model.solve().f * 4) / (-optswap.model.reactions.EX_glc__D_e.flux * 6)
with TimeMachine() as tm:
    optswap.model.swap_reaction('GAPD', tm)
    swapped = (optswap.model.solve().f * 4) / (-optswap.model.reactions.EX_glc__D_e.flux * 6)

print("product/substrate yield without swap: {}\nproduct/substrate yield with swap: {}".format(original, swapped))