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

# # Simulating with FBA
# 
# Simulations using flux balance analysis can be solved using Model.optimize(). This will maximize or minimize (maximizing is the default) flux through the objective reactions.

# In[1]:


import pandas
pandas.options.display.max_rows = 100

import cobra.test
model = cobra.test.create_test_model("textbook")


# ## Running FBA

# In[2]:


model.optimize()


# The Model.optimize() function will return a Solution object, which will also be stored at model.solution. A solution object has several attributes:
# 
#  - f: the objective value
#  - status: the status from the linear programming solver
#  - x_dict: a dictionary of {reaction_id: flux_value} (also called "primal")
#  - x: a list for x_dict
#  - y_dict: a dictionary of {metabolite_id: dual_value}.
#  - y: a list for y_dict

# For example, after the last call to model.optimize(), the status should be 'optimal' if the solver returned no errors, and f should be the objective value

# In[3]:


model.solution.status


# In[4]:


model.solution.f


# ### Analyzing FBA solutions
# Models solved using FBA can be further analyzed by using summary methods, which output printed text to give a quick representation of model behavior. Calling the summary method on the entire model displays information on the input and output behavior of the model, along with the optimized objective.

# In[5]:


model.summary()


# In addition, the input-output behavior of individual metabolites can also be inspected using summary methods. For instance, the following commands can be used to examine the overall redox balance of the model

# In[6]:


model.metabolites.nadh_c.summary()


# Or to get a sense of the main energy production and consumption reactions

# In[7]:


model.metabolites.atp_c.summary()


# ## Changing the Objectives
# 
# The objective function is determined from the objective_coefficient attribute of the objective reaction(s). Generally, a "biomass" function which describes the composition of metabolites which make up a cell is used.

# In[8]:


biomass_rxn = model.reactions.get_by_id("Biomass_Ecoli_core")


# Currently in the model, there is only one objective reaction (the biomass reaction), with an objective coefficient of 1.

# In[9]:


model.objective


# The objective function can be changed by assigning Model.objective, which can be a reaction object (or just it's name), or a dict of {Reaction: objective_coefficient}.

# In[10]:


# change the objective to ATPM
model.objective = "ATPM"

# The upper bound should be 1000, so that we get
# the actual optimal value
model.reactions.get_by_id("ATPM").upper_bound = 1000.
model.objective


# In[11]:


model.optimize().f


# The objective function can also be changed by setting Reaction.objective_coefficient directly.

# In[12]:


model.reactions.get_by_id("ATPM").objective_coefficient = 0.
biomass_rxn.objective_coefficient = 1.

model.objective


# ## Running FVA
# 
# FBA will not give always give unique solution, because multiple flux states can achieve the same optimum. FVA (or flux variability analysis) finds the ranges of each metabolic flux at the optimum.

# In[13]:


fva_result = cobra.flux_analysis.flux_variability_analysis(
    model, model.reactions[:20])
pandas.DataFrame.from_dict(fva_result).T.round(5)


# Setting parameter fraction_of_optimium=0.90 would give the flux ranges for reactions at 90% optimality.

# In[14]:


fva_result = cobra.flux_analysis.flux_variability_analysis(
    model, model.reactions[:20], fraction_of_optimum=0.9)
pandas.DataFrame.from_dict(fva_result).T.round(5)


# ### Running FVA in summary methods
# Flux variability analysis can also be embedded in calls to summary methods. For instance, the expected variability in substrate consumption and product formation can be quickly found by

# In[15]:


model.optimize()
model.summary(fva=0.95)


# Similarly, variability in metabolite mass balances can also be checked with flux variability analysis

# In[16]:


model.metabolites.pyr_c.summary(fva=0.95)


# In these summary methods, the values are reported as a the center point +/- the range of the FVA solution, calculated from the maximum and minimum values.

# ## Running pFBA
# 
# Parsimonious FBA (often written pFBA) finds a flux distribution which gives the optimal growth rate, but minimizes the total sum of flux. This involves solving two sequential linear programs, but is handled transparently by cobrapy. For more details on pFBA, please see [Lewis et al. (2010)](http://dx.doi.org/10.1038/msb.2010.47).

# In[17]:


FBA_sol = model.optimize()
pFBA_sol = cobra.flux_analysis.optimize_minimal_flux(model)


# These functions should give approximately the same objective value

# In[18]:


abs(FBA_sol.f - pFBA_sol.f)