Stages¶

import impact
import cobra
import cobra.test

# Hide some warnings to make the output cleaner
import warnings
warnings.filterwarnings('ignore')

% matplotlib inline

# Let's grab the iJO1366 E. coli model, cobra's test module has a copy
model = cobra.test.create_test_model("ecoli")

# Optimize the model
model.optimize()

# Print a summary of the fluxes
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
# We'll use numpy to generate an arbitrary time vector
import numpy as np

# The initial conditions (mM) [biomass, substrate,
#                              product1, product2, ..., product_n]
y0 = [0.05, 1000, 0, 0, 0, 0]
t_aerobic = np.linspace(0,6,600)

# Returns a dictionary of the profiles
from impact.synthetic_data import generate_data
dFBA_profiles_aerobic = generate_data(y0, t_aerobic, model,
                                      biomass_keys, substrate_keys,
                                      product_keys, plot = True)

# Let's make a new vector of initial values based on the output of the previous simulation
y0 = [dFBA_profiles_aerobic[analyte][-1] for analyte in biomass_keys+substrate_keys+product_keys]

# Optimize the model after simulating anaerobic conditions,
# we don't get many products aerobically
model.reactions.get_by_id('EX_o2_e').lower_bound = 0
model.reactions.get_by_id('EX_o2_e').upper_bound = 0
model.optimize()
model.summary()

t_anaerobic = np.linspace(6,10,400)
dFBA_profiles_anaerobic = generate_data(y0, t_anaerobic, model,
                                        biomass_keys, substrate_keys,
                                        product_keys, plot = True)

# Plot the whole profile
combined_profiles = {analyte:
                     dFBA_profiles_aerobic[analyte]+
                     dFBA_profiles_anaerobic[analyte]
                     for analyte in biomass_keys+substrate_keys+product_keys}

import matplotlib.pyplot as plt
% matplotlib inline
plt.figure(figsize=[12,6])
for analyte in substrate_keys + product_keys:
    plt.plot(np.linspace(0,10,1000),combined_profiles[analyte])
    min,max = plt.ylim()
    plt.ylim([0,max])

plt.text(5,max*0.95,'Aerobic')
plt.text(6+.5,max*0.95,'Anaerobic')
plt.legend(substrate_keys+product_keys,loc=2)
plt.plot([6]*10,np.linspace(min,max,10),'k',linewidth=3)

plt.figure(figsize=[12,6])
plt.plot(np.linspace(0,10,1000),combined_profiles[biomass_keys[0]])
min,max = plt.ylim()
plt.ylim([0,max])
plt.plot([6]*10,np.linspace(min,max,10),'k',linewidth=3)
plt.text(5,max*0.95,'Aerobic')
plt.text(6+.5,max*0.95,'Anaerobic')
plt.legend(['Biomass'],loc=2)

# These time courses together form a single trial
single_trial = impact.SingleTrial()
for analyte in analyte_keys:
    # Instantiate the timecourse
    timecourse = impact.TimeCourse()

    # Define the trial identifier for the experiment
    timecourse.trial_identifier.analyte_name = analyte
    timecourse.trial_identifier.strain_id = 'Demo'

    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 = np.linspace(0,10,1000)
    timecourse.data_vector = combined_profiles[analyte]
    single_trial.add_titer(timecourse)

# 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)

import matplotlib.pyplot as plt
for single_trial in replicate_trial.single_trial_list:
    print(replicate_trial.single_trial_list[0].analyte_dict[biomass_keys[0]].rate)
    plt.plot(np.linspace(0,10,1000),combined_profiles[biomass_keys[0]])
    plt.plot([0,1,2,3,4,5,6,7,8,9,10],replicate_trial.single_trial_list[0].analyte_dict[biomass_keys[0]].data_curve_fit([0,1,2,3,4,5,6,7,8,9,10]))

# We can break the trials into stages
replicate_trial.calculate_stages(stage_indices=[[0,600],[600,1000]])

for stage_number in [0,1]:
    print('Stage: ',stage_number)
    plt.figure()
    stage = replicate_trial.stages[stage_number]
    timecourse = stage.single_trial_list[0].analyte_dict[biomass_keys[0]]

    for single_trial in replicate_trial.single_trial_list:
        print(timecourse.rate)
        plt.plot(timecourse.time_vector,
                timecourse.data_vector)
        plt.plot(timecourse.time_vector,timecourse.data_curve_fit(timecourse.time_vector))

    plt.show()