<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Analyzing Features
Rates
Step1:  Now that we have a simulated 'two-stage' fermentation data, let's try to analyze it. First let's try to curve fit and pull parameters from the overall data
Step2:  Specific Productivity
The specific productivity is the rate, normalized to biomass concentration. This is the value used by dFBA to calculate the initial growth curves.
$$ \frac{dX}{dt} = \mu X $$
$$ \frac{dP}{dt} = \dot{P} X $$
Here we find the median specific productivity and compare it to the FBA solutions
Step3:  Here we can see that the specific productivities match those calculated in the COBRA model (which is expected since that's where the generated data came from). For cells in steady-state, this should be a constant.
Metabolic Model Integration
With specific productivities, metabolic model integration becomes straightforward. We can constrain the export fluxes based on the data, and solve the missing fluxes. Since this is model-generated data, this will close the mass balance entirely.
For experimental data, this offers a prediction to help close the balance.
Step4:  dFBA Batch Integration
We can integrate the experimental data across the course of the batch, and solve an FBA for each point to account for any misisng metabolites and get a dynamic flux map for the batch.
Step5:  Carbon Balance
From here, the carbon balance is straight forward. The export fluxes match experimental data, the remaining fluxes are estimated, and we can complete the carbon balance. Units are import here, and should all be converted a standard mass unit (g/L) to complete the balance.
Step6:  Flux mapping with escher
More detail is provided in the escher documentation
<ASSISTANT_TASK:>
Python Code:
import impact as impt
import cobra
import cobra.test
import cobra.io
import numpy as np
# import matplotlib.pyplot as plt
% matplotlib inline
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
from plotly.graph_objs import Bar, Layout, Figure, Scatter
init_notebook_mode()
# We include this to ensure the js is loaded when viewed online
# from IPython.display import HTML
# HTML('<script src="https://cdn.plot.ly/plotly-latest.min.js"></script>')
import impact.plotting as implot
# Let's grab the iJO1366 E. coli model, cobra's test module has a copy
model = cobra.test.create_test_model("ecoli")
# Simulate anaerobic conditions by prevent oxygen uptake
model.reactions.get_by_id('EX_o2_e').knock_out() 
# Optimize the model
sol = model.optimize()
model.summary()
# Let's consider one substrate and five products
biomass_keys = ['Ec_biomass_iJO1366_core_53p95M']
substrate_keys = ['EX_glc_e']
product_keys = ['EX_for_e','EX_ac_e','EX_etoh_e','EX_succ_e']
analyte_keys = biomass_keys+substrate_keys+product_keys
# The initial conditions (mM) [biomass, substrate, 
#                              product1, product2, ..., product_n]
y0 = [0.05, 100, 0, 0, 0, 0]
t = np.linspace(0,16,8)
# Returns a dictionary of the profiles
from impact.helpers.synthetic_data import generate_data
dFBA_profiles = generate_data(y0, t, model, 
                                      biomass_keys, substrate_keys, 
                                      product_keys, plot = False)
implot.plot([implot.go.Scatter(x=t,y=dFBA_profiles[exchange],name=exchange) for exchange in dFBA_profiles])
Explanation: Analyzing Features
Rates
End of explanation
# These time courses together form a single trial
single_trial = impt.SingleTrial()
for analyte in analyte_keys:
    # Instantiate the timecourse
    timecourse = impt.TimeCourse()
    
    # Define the trial identifier for the experiment
    ti = impt.ReplicateTrialIdentifier()
    ti.strain.nickname = 'my_strain'
    ti.analyte_name = analyte
    timecourse.trial_identifier = ti
    
    if analyte in biomass_keys:
        timecourse.trial_identifier.analyte_type = 'biomass'
    elif analyte in substrate_keys:
        timecourse.trial_identifier.analyte_type = 'substrate'
    elif analyte in product_keys:
        timecourse.trial_identifier.analyte_type = 'product'
    else:
        raise Exception('unid analyte')
    
    timecourse.time_vector = t
    timecourse.data_vector = dFBA_profiles[analyte]
    single_trial.add_analyte_data(timecourse)
single_trial.calculate()
# # Add this to a replicate trial (even though there's one replicate)
# replicate_trial = impact.ReplicateTrial()
# replicate_trial.add_replicate(single_trial)
# # Add this to the experiment
# experiment = impact.Experiment(info = {'experiment_title' : 'test experiment'})
# experiment.add_replicate_trial(replicate_trial)
# import plotly.offline
# plotly.offline.init_notebook_mode()
# fileName = impact.printGenericTimeCourse_plotly(replicateTrialList=[replicate_trial], 
#                                                 titersToPlot=biomass_keys,` output_type='image',)
# from IPython.display import Image
# Image(fileName)
# We can query the rates instead.
for keys in [biomass_keys,substrate_keys,product_keys]:
#     plt.figure()
    data = []
    for analyte in keys:
        trace = Scatter(x=t,
            y=single_trial.analyte_dict[analyte].gradient,
            name = analyte)
        data.append(trace)
    fig = Figure(data=data,
                 layout=Layout(title="Calculated rates", 
                               yaxis = {'title':'Production rate [mmol /  h]'})
                 )
    iplot(fig)
#         plt.plot(t,single_trial.analyte_dict[analyte].gradient)
#     plt.legend(keys)
Explanation: Now that we have a simulated 'two-stage' fermentation data, let's try to analyze it. First let's try to curve fit and pull parameters from the overall data
End of explanation
exchange_keys = biomass_keys + substrate_keys + product_keys
median_specific_productivity = {}
model_specific_productivity = {}
for analyte in exchange_keys:
    model_specific_productivity[analyte] = model.solution.x_dict[analyte]
    median_specific_productivity[analyte] = np.median(single_trial.analyte_dict[analyte].specific_productivity.data)
trace = Bar(x=[analyte for analyte in median_specific_productivity],
            y=[median_specific_productivity[analyte] for analyte in median_specific_productivity],
           name = 'Calculated')
trace2 = Bar(x=[analyte for analyte in model_specific_productivity],
            y=[model_specific_productivity[analyte] for analyte in model_specific_productivity],
            name = 'Model')
data = [trace, trace2]
fig = Figure(data=data,
             layout=Layout(title="Calculated vs Model Specific Productivity", 
                           yaxis = {'title':'Specific Productivity [mmol / (gdw h)]'})
             )
iplot(fig)
Explanation: Specific Productivity
The specific productivity is the rate, normalized to biomass concentration. This is the value used by dFBA to calculate the initial growth curves.
$$ \frac{dX}{dt} = \mu X $$
$$ \frac{dP}{dt} = \dot{P} X $$
Here we find the median specific productivity and compare it to the FBA solutions
End of explanation
# Let's take the model and add bounds for the known reactions
for analyte in analyte_keys:
    with_noise = np.median(single_trial.analyte_dict[analyte].specific_productivity.data)
    print(analyte,' ',with_noise)
    model.reactions.get_by_id(analyte).lower_bound = with_noise
    model.reactions.get_by_id(analyte).upper_bound = with_noise
print('\n\n\n')
solution = model.optimize()
model.summary()
Explanation: Here we can see that the specific productivities match those calculated in the COBRA model (which is expected since that's where the generated data came from). For cells in steady-state, this should be a constant.
Metabolic Model Integration
With specific productivities, metabolic model integration becomes straightforward. We can constrain the export fluxes based on the data, and solve the missing fluxes. Since this is model-generated data, this will close the mass balance entirely.
For experimental data, this offers a prediction to help close the balance.
End of explanation
# dFBA based batch simulation
# import copy
# model2 = copy.deepcopy(model)
# impact.helpers.synthetic_data.dynamic_model_integration(t, y0, model2, single_trial, biomass_keys, substrate_keys, product_keys, extra_points_multiplier = 5)
Explanation: dFBA Batch Integration
We can integrate the experimental data across the course of the batch, and solve an FBA for each point to account for any misisng metabolites and get a dynamic flux map for the batch.
End of explanation
# Convert mmol/gdw - hr to g / gdw - hr
# Create a dictionary for the mapping of flux reactions to metabolite names to get the molar mass
metabolite_export_reaction_dict = {
    'EX_for_e':'for_e',
    'EX_ac_e':'ac_e',
    'EX_etoh_e':'etoh_e',
    'EX_glc_e':'glc__D_e',
    'EX_succ_e':'succ_e'
}
mass_flux = {}
for analyte in product_keys+substrate_keys:
    mass_flux[analyte] = model.solution.x_dict[analyte] \
    * model.metabolites.get_by_id(metabolite_export_reaction_dict[analyte]).formula_weight
mass_flux[biomass_keys[0]] = model.solution.x_dict[biomass_keys[0]]*1000
print(mass_flux)
# The balance is the sum of all metabolites, the uptake is already negative.
balance = sum(mass_flux[metabolite] for metabolite in mass_flux)
print(balance)
# The closure is the total substrate accounted for
percent_closure = balance/mass_flux['EX_glc_e']
print('\nThe mass balance is %f%% closed' % ((1-percent_closure)*100))
print(mass_flux)
# Two situations
# 1: balance is > 100% (must be some unaccounted substrate)
labels = [analyte.split('_')[1] for analyte in mass_flux]
print(percent_closure)
if percent_closure < 0:
    fig = {
    'data': [{'labels': labels + ['Unaccounted substrate'],
              'values': [mass_flux[metabolite] for metabolite in mass_flux]+[balance],
              'type': 'pie'}],
    'layout': {'title': 'Mass balance'}
     }
# 2: balance is < 100% (must be some unaccounted product)
else:
    fig = {
    'data': [{'labels': labels + ['Missing product'],
              'values': [mass_flux[metabolite] for metabolite in mass_flux]+[abs(balance)],
              'type': 'pie'}],
    'layout': {'title': 'Mass balance'}
         }
iplot(fig)
analyte_keys = substrate_keys + biomass_keys + product_keys
# Calculate the waterfall for the products and biomass
substrate_consuming_keys = biomass_keys+product_keys
running_total = balance
base = []
for i, key in enumerate(analyte_keys):
    if i == 0:
        running_total += -mass_flux[key]
        base.append(0)
    else:
        base.append(running_total)
        running_total += -mass_flux[key]
    
base_trace = Bar(x=analyte_keys,
                 y=base,
                 marker=dict(color='rgba(1,1,1, 0.0)',))
substrate_trace = Bar(x=[substrate_keys[0]],
                      y=[-mass_flux[substrate_keys[0]]],
                        marker=dict(color='rgba(55, 128, 191, 0.7)',
                                    line=dict(color='rgba(55, 128, 191, 1.0)',
                                              width=2,)))
biomass_trace = Bar(x=[biomass_keys[0]],
                    y=[-mass_flux[biomass_keys[0]]],
                       marker=dict(color='rgba(50, 171, 96, 0.7)',
                                   line=dict(color='rgba(50, 171, 96, 1.0)',
                                             width=2,)))
products_trace = Bar(x=product_keys,
                     y=[-mass_flux[metabolite] for metabolite in product_keys],
                      marker=dict(color='rgba(219, 64, 82, 0.7)',
                                  line=dict(color='rgba(219, 64, 82, 1.0)',
                                            width=2,)),
                    )
balance_base_trace = Bar(x=['Balance'],
                 y=[base[1]-balance],
                 marker=dict(color='rgba(1,1,1, 0.0)',))
balance_trace = Bar(x=['Balance'],
                     y=[balance],
                      marker=dict(color='rgba(0, 0, 0, 0.7)',
                                  line=dict(color='rgba(0, 0, 0, 1.0)',
                                            width=2,)),
                    )
data = [substrate_trace, balance_base_trace, base_trace, balance_trace, biomass_trace,products_trace]
layout = Layout(barmode='stack',
#                 paper_bgcolor='rgba(245, 246, 249, 1)',
#                 plot_bgcolor='rgba(245, 246, 249, 1)',
                yaxis={'title': 'Exchange rate (mmol/(gdw h))'},
                showlegend=False
            )
fig = Figure(data=data,layout=layout)
iplot(fig)
Explanation: Carbon Balance
From here, the carbon balance is straight forward. The export fluxes match experimental data, the remaining fluxes are estimated, and we can complete the carbon balance. Units are import here, and should all be converted a standard mass unit (g/L) to complete the balance.
End of explanation
import escher
import json
from IPython.display import HTML
b = escher.Builder(map_name='iJO1366.Central metabolism',
                   reaction_data=solution.x_dict,
                   # color and size according to the absolute value
                   reaction_styles=['color', 'size', 'abs', 'text'],
                   # change the default colors
                   reaction_scale=[{'type': 'min', 'color': '#cccccc', 'size': 4},
                                   {'type': 'mean', 'color': '#0000dd', 'size': 20},
                                   {'type': 'max', 'color': '#ff0000', 'size': 40}],
                   # only show the primary metabolites
                   hide_secondary_metabolites=True)
b.display_in_notebook()
Explanation: Flux mapping with escher
More detail is provided in the escher documentation: http://nbviewer.jupyter.org/github/zakandrewking/escher/blob/master/docs/notebooks/COBRApy%20and%20Escher.ipynb
End of explanation
<END_TASK>