<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Tellurium Introduction
Step1:  Everything in a single tool - tellurium
Step2:  Antimony is a language that is analog to SBML Systems Biology Markup Language but human-readable. Therefore the concepts that are present in SBML can be found in antimony too. The definition of a compartment is such a SB concept!
Step 1
Step3:  Ok... But that's just boring. Let's define something truly special and add a second compartment which is located in your previous one
Step4:  Awesome! But be aware... circular definitions are not allowed.
Something like this won't work
Step5:  To specify units you can write the value and the unit combined in the same line
Step6:  Or simply change the units
Step7:  Step 2
Step8:  Smart! You already initialized your species. But you forgot to specify the location of your species! Add these lines
Step9:  Step 3
Step10:  Syntax in antimony is subtile
Step11:  Or set a time dependend event. You may notice that time is special and the already predefined model time
Step12:  Ok, now you may want to test your model and simulate it. You should definitly do this before we continue. Make sure that your model is correctly defined at the top of this jupyter notebook (you may execute the cell again!).
Then you can run it
Step13:  Import/Export Models
Step14:  Load a model... You can choose from
Step15:  ii) download 'curated BioModels translated to Antimony' from here http
Step16:  iii) load your own model
Step17:  You can also convert models between modeling-languages
Step18:  Then save the SBML model in a file
Step19:  An read it again...
Step20:  You can also draw your model
Step21:  Parameterscan
Important concept
Step22:  Stochastic vs. deterministic
<ASSISTANT_TASK:>
Python Code:
import tellurium as te; te.setDefaultPlottingEngine('matplotlib')
%matplotlib inline
antimony_model = '''J0: -> y; -x;J1: -> x; y;x = 1.0;y = 0.2;'''
r = te.loada(antimony_model)
r.simulate(0,100,1000)
r.plot()
Explanation: Tellurium Introduction:
Motivation ... a minimal example!
Just a few lines of code allow a model to be loaded, simulated and plotted!
End of explanation
import tellurium as te
model = ''''''
model_backup = '''
model example
    # UNITS:
    #
    #unit alienliters = 0.123 liters
    #unit fL = 1.e-15 liters
    # COMPARTMENTS:
    #
    compartment cell = 100;
    compartment mitochondria = 10;
    mitochondria in cell;
    #mitochondria has fL;
    #cell has fL;
    # INITIAL VALUES
    #
    TOM1 = 10;
    A_c = 100;
    A_m = 1;
    I = 1;
    TOM1 in cell;
    A_c in cell;
    A_m in mitochondria;
    I in cell;
    # REACTIONS
    #
    T0: A_c + TOM1 -> 2 A_m + TOM1; kineticLaw;
    kineticLaw := k1 * TOM1 * (A_c - A_m)/I
    k1 = 0.01
    # EVENTS:
    #E1: at (A_c < 50): k1 = 0;
    #E2: at (time>300): TOM1 = 20, A_c = 120;
end
'''
r = te.loada(model)
def aa(line, cell):
    global model
    #print line,cell
    model = model + str(cell) + "\n"
    te.loada(str(model))
get_ipython().register_magic_function(aa, "cell")     
Explanation: Everything in a single tool - tellurium:
As you could see, tellurium allows you to create and manipulate a model in a pythonic way. The tellurium package consists of formerly two distinct packages: antimony which is a modeling-language and roadrunner a sbml-simulator. The main advantage of tellurium is how both work together in python: as soon as you load a model a roadrunner simulator will be instanciated simultaniously. This way your model and the solver are a python object - the same object!  This allows you to have real power over your model and the solver and control everything within your python script. 
For further reading:
Full documentation
http://tellurium.readthedocs.io/en/latest/index.html
Short documentation
http://tellurium.analogmachine.org/documentation/antimony-documentation
Antimony (bio-)models library
http://antimony.sourceforge.net/antimony-biomodels.html
Tutorial Part I - Create a Model :
Let us start from scratch and define a new model (if you execute the next cell then tellurium tries to load your model and performs a basic syntax check at the same time!):
End of explanation
%%aa # %%aa: A-ppend to A-ntimony model, only for educational purpose!
compartment cell = 100;
Explanation: Antimony is a language that is analog to SBML Systems Biology Markup Language but human-readable. Therefore the concepts that are present in SBML can be found in antimony too. The definition of a compartment is such a SB concept!
Step 1: Define a compartment!
End of explanation
%%aa
compartment mitochondria;
mitochondria = 10;
mitochondria in cell;
Explanation: Ok... But that's just boring. Let's define something truly special and add a second compartment which is located in your previous one:
End of explanation
%%aa
unit alienliters = 0.123 liters
unit fL = 1.e-15 liters
Explanation: Awesome! But be aware... circular definitions are not allowed.
Something like this won't work:
comp_a in comp_b; comp_b in comp_c; comp_c in comp_a;
SBML has a system of unit definition. But the numbers we assigned to the compartments are without any specified unit at the moment. It is good practice to always use units and it is possible to define custom units!
End of explanation
%%aa
mitochondria = 10 fL;
Explanation: To specify units you can write the value and the unit combined in the same line:
End of explanation
%%aa
mitochondria has fL;
Explanation: Or simply change the units:
End of explanation
%%aa
TOM1 = 10;
A_c = 100;
A_m = 1;
Explanation: Step 2: Define some species.
End of explanation
%%aa
TOM1 in cell;
A_c in cell;
A_m in mitochondria;
Explanation: Smart! You already initialized your species. But you forgot to specify the location of your species! Add these lines:
End of explanation
%%aa
T0: A_c + TOM1 -> A_m + TOM1; kineticLaw;
kineticLaw := k1 * TOM1 * (A_c - (A_m+0.01))
k1 = 0.01
Explanation: Step 3: Define reactions.
End of explanation
%%aa
E1: at (A_c < 50): k1 = 0;
Explanation: Syntax in antimony is subtile:
There are 3 (+1) important bits in a reaction:
0) The name of your reaction:
name: A + E -> B + 2C + E; kineticLaw;
i) The direction '->' (sets the sign in an ODE in front the kineticLaw: plus or minus):
name: A + E -> B + 2C + E; kineticLaw;
ii) The stoichiometry:
name: A + E -> B + 2C + E; kineticLaw;
iii) The kinetic law:
name: A + E -> B + 2C + E; kineticLaw;
The kinetic law is always in units of amount/time! This way it is consistent with stochastic simulations (you can directly run the 'gillespie'-algorithm in tellurium's solver: roadrunner) and furthermore any considerations regarding compartment sizes (even time depending compartment changes or transport reactions) are automatically done by the solver to get the correct concentrations for your species.
About Assignments and ODE's:
This " = " is a simple assignment, basically an initialization:
parameter = 5
species = 1
species2 = species + parameter
This " := " is a time-dependend assignment, the variable at the left will adapt, if the right hand side changes:
species0 := 3 * time
species1 := species0
This " '= " will create an ODE:
species'= k * species
Step 4: Events.
You can define events, which may change values depending on boolean expressions...
End of explanation
%%aa
E2: at (time>300): TOM1 = 20, A_c = 120;
Explanation: Or set a time dependend event. You may notice that time is special and the already predefined model time:
End of explanation
import tellurium as te
r = te.loada(model)
r.simulate(0,1000,1000)
r.plot()
Explanation: Ok, now you may want to test your model and simulate it. You should definitly do this before we continue. Make sure that your model is correctly defined at the top of this jupyter notebook (you may execute the cell again!).
Then you can run it:
End of explanation
import tellurium as te
Explanation: Import/Export Models:
First, we need the tellurium package:
End of explanation
print te.listTestModels()
#r_feedback = te.loadTestModel('feedback.xml')
#result = r_feedback.simulate()
#r_feedback.plot()
Explanation: Load a model... You can choose from:
i) load a test model from tellurium itself.
End of explanation
import urllib2
bio1 = urllib2.urlopen('http://antimony.sourceforge.net/examples/biomodels/BIOMD0000000001.txt').read()
r_bio1 = te.loadAntimonyModel(bio1)
r_bio1.simulate(0,1000,10000)
r_bio1.plot()
Explanation: ii) download 'curated BioModels translated to Antimony' from here http://antimony.sourceforge.net/antimony-biomodels.html:
End of explanation
import tellurium as te
r = te.loada(model)
r.simulate()
r.plot()
Explanation: iii) load your own model:
End of explanation
sbml_model = te.antimonyToSBML(model)
Explanation: You can also convert models between modeling-languages:
End of explanation
with open('model.xml','wb') as f:
	f.write(sbml_model)
# directly for antimony:
with open('example.antimony','wb') as f:
    f.write(model)
	#f.write(model_backup)
Explanation: Then save the SBML model in a file:
End of explanation
sbml_model = te.loadSBMLModel('model.xml')
example = te.loada('example.antimony')
Explanation: An read it again...
End of explanation
def draw(model):
    diagram = te.visualization.SBMLDiagram(model.getSBML())
    diagram.draw()
draw(r)
Explanation: You can also draw your model:
End of explanation
import numpy as np
import matplotlib.pyplot as plt
r = te.loada(model)
fig = plt.figure(figsize=(10,10))
for i in np.linspace(0.001,0.005,5):
    r.reset()
    r.k1 = i    
    #r.reset() # reset species to initial values
    #r.resetAll() # reset species + paramters
    #r.resetToOrigin() # similar to load the model
    r.integrator='cvode'
    result = r.simulate(0,1000)  
    plt.plot(result[:,0],result[:,1:],label='k1=%s'%i)
    
plt.legend()
plt.show()
Explanation: Parameterscan
Important concept: Reset a model!
.reset() # reset species to initial values 
.resetAll() # reset species + paramters
.resetToOrigin() # similar to load the model
End of explanation
r = te.loada(model_backup)
    
#plt.close()
fig = plt.figure(figsize=(10,10))
r.resetToOrigin()
r.timeCourseSelections=['time','[A_c]']
for i in range(10):
    r.reset()    
    r.integrator = 'gillespie'
    result = r.simulate(0,100,1000)
    plt.plot(result[:,0],result[:,1:],color = 'lightblue',label='gillespie sim')
    
r.reset()
r.integrator='cvode'
result = r.simulate(0,100,1000)   
plt.plot(result[:,0],result[:,1:],color = 'red',label='cvode sim')
      
plt.legend()
plt.show()
r = te.loada(model)
#print r.getAntimony()
model = '''
// Created by libAntimony v2.9.4
// Compartments and Species:
compartment cell, mitochondria;
species TOM1 in cell, A_c in cell, A_m in mitochondria;
// Assignment Rules:
kineticLaw := k1*TOM1*(A_c - (A_m + 0.01));
// Reactions:
T0: A_c + TOM1 -> A_m + TOM1; kineticLaw;
// Events:
E1: at A_c < 50: k1 = 0;
E2: at time > 300: TOM1 = 20, A_c = 120;
// Species initializations:
TOM1 = 10;
A_c = 100;
A_m = 1;
// Compartment initializations:
cell = 100;
mitochondria = 10;
mitochondria has fL;
// Variable initializations:
k1 = 0.01;
// Other declarations:
var k1, kineticLaw;
const cell, mitochondria;
// Unit definitions:
unit alienliter = 1.23e-1 litre;
unit fL = 1e-15 litre;
const A_c; 
'''
r = te.loada(model)
r.timeCourseSelections= ['time'] + ['['+ r.getBoundarySpeciesIds()[0]+ ']'] + r.getFloatingSpeciesIds()
r.simulate(0,1000)
r.plot()
Explanation: Stochastic vs. deterministic
End of explanation
<END_TASK>