In [ ]:

    

%matplotlib notebook
 
from libsignetsim import SbmlDocument, Experiment, TimeseriesSimulation, ModelVsTimeseriesOptimization

from os.path import join
import matplotlib.pyplot as plt


    

In [ ]:

    

# Loading the model
sbml_filename = join("sbml_files", "ras_mapk.xml")
sbml_doc = SbmlDocument()
sbml_doc.readSbmlFromFile(sbml_filename)
model = sbml_doc.getModelInstance()

# Putting FGF2 stimulus as initial condition, otherwise it's boring
model.listOfSpecies.getBySbmlId("fgf2").setValue(333)


    

In [ ]:

    

# Simulating the model as is
sim = TimeseriesSimulation([model], time_min=0, time_ech=6, time_max=3600)
sim.run()
ts, ys = sim.getRawData()[0]


    

In [ ]:

    

# Generating the fake experimental data from the simulation
experiment = Experiment()
condition = experiment.createCondition()

for i, t in enumerate(ts):
    condition.addObservation(t, "Ras-GTP", ys["ras_gtp"][i])
    condition.addObservation(t, "Mapk-PP", ys["mapk_pp"][i])


    

In [ ]:

    

# Modifying one parameter
param = model.listOfParameters.getBySbmlId("raf_activation_cat")
print "Original parameter value : %g" % param.getValue()
param.setValue(1)


    

In [ ]:

    

# Simulating the model with this modification
sim_mod = TimeseriesSimulation([model], time_min=0, time_ech=6, time_max=3600)
sim_mod.run()
ts_mod, ys_mod = sim_mod.getRawData()[0]


    

In [ ]:

    

# Comparing the two trajectories for Ras-GTP and MAPK-PP
plt.figure()
plot1 = plt.subplot(2,1,1)
plot1.plot(ts, ys["ras_gtp"], ts, ys["mapk_pp"])
plot2 = plt.subplot(2,1,2)
plot2.plot(ts_mod, ys_mod["ras_gtp"], ts_mod, ys_mod["mapk_pp"])


    

In [ ]:

    

# Selecting which parameter to adjust in the optization, its initial value, and its boundaries
selected_parameters = []
selected_parameters.append((param, 1, 1e-8, 1e+8))

# Doing the fitting
# Lambda is a quality parameter. The smaller, the longer, but the highest chance to get the global optimum.
# Depends on the complexity of the problem.

# Freeze count is another quality parameter, only used here to make a quick example, since the problem is really easy.
# It defined how many optimization iterations to go through before stopping, if the score doesn't evolve anymore. 
# Default is 100. 
fit = ModelVsTimeseriesOptimization(
    workingModel=model,
    list_of_experiments=[experiment],
    parameters_to_fit=selected_parameters,
    p_lambda=0.01, p_freeze_count=10
)

score = fit.runOptimization(2)
parameters = fit.readOptimizationOutput()


    

In [ ]:

    

# Settings the parameter to the fitted value, and simulating again
fitted_value = parameters[param]
print "> Fitted parameter value = %g" % fitted_value
param.setValue(fitted_value)
sim_fitted = TimeseriesSimulation([model], time_min=0, time_ech=6, time_max=3600)
sim_fitted.run()
ts_fitted, ys_fitted = sim_fitted.getRawData()[0]


    

In [ ]:

    

# Comparing the fitting result to the original.
fig2 = plt.figure()
plot1 = plt.subplot(2,1,1)
plot1.plot(ts, ys["ras_gtp"], ts, ys["mapk_pp"])
plot2 = plt.subplot(2,1,2)
plot2.plot(ts_fitted, ys_fitted["ras_gtp"], ts_fitted, ys_fitted["mapk_pp"])