# 3- ODE solving (dynamics.py module).¶

This chapter shows how to use 4 of the most common S-timator functions:

• read_model(), reads a string that conforms to a model description language, returning a Model object
• solve(), computes a solution of the ODE system associated with a model.
• scan(), calls Model.solve() several times, scanning a model parameter in a range of values.
• plot(), draws a graph of the results returned from solve() or scan().
from stimator import read_model, examples, __version__
print "S-timator version", __version__.fullversion, "(%s)" %__version__.date

S-timator version 0.9.95 (Jan 2015)


## Example 1¶

Glyoxalase system model

mdl = examples.models.glyoxalases.text
print mdl

s = m1.solve(tf=4030.0)
s.plot()

print '==== Last time point ===='
print 'At t =', s.t[-1]
for x in s.last:
print "%-8s= %f" % (x, s.last[x])

title Glyoxalase system
glo1 = HTA -> SDLTSH  , rate = V1*HTA/(Km1 + HTA)
glo2 = SDLTSH -> DLac , rate = V2*SDLTSH/(Km2 + SDLTSH)
V1  = 2.57594e-05
Km1 = 0.252531
V2  = 2.23416e-05
Km2 = 0.0980973
init :(SDLTSH = 7.69231E-05, HTA = 0.1357)
tf: 4030

==== Last time point ====
At t = 4030.0
SDLTSH  = 0.022744
DLac    = 0.010462
HTA     = 0.102571


## Example 2¶

Branched pathway

from numpy import append, linspace
mdl = examples.models.branched.text; print mdl

times = append(linspace(0.0, 5.0, 500), linspace(5.0, 10.0, 500))

m2.solve(tf=10.0, times=times).plot()

title Branched pathway

v1 = A -> B, k1*A,      k1 = 10
v2 = B -> C, k2*B**0.5, k2 = 5
v3 = C -> D, k3*C**0.5, k3 = 2
v4 = C -> E, k4*C**0.5, k4 = 8
v5 = D ->  , k5*D**0.5, k5 = 1.25
v6 = E ->  , k6*E**0.5, k6 = 5
A  = 0.5
init : (B = 2, C = 0.25, D = 0.64, E = 0.64)
!! B D
tf: 10


## Example 3¶

Calcium spikes: CICR model

mdl = examples.models.ca.text; print mdl
#chaining functions...

title CICR model. Calcium spikes
v0         = -> Ca, 1
v1         = -> Ca, k1*B*step(t, t_stimulus)

k1         = 7.3
B          = 0.4
t_stimulus = 1.0

export     = Ca ->  , 10 ..
leak       = CaComp -> Ca, 1 ..

v2         = Ca -> CaComp, 65 * Ca**2 / (1+Ca**2)
v3         = CaComp -> Ca, 500*CaComp**2/(CaComp**2+4) * Ca**4/(Ca**4 + 0.6561)
init       : Ca = 0.1, CaComp = 0.63655


## Example 4¶

Rossler chaotic system

mdl = examples.models.rossler.text; print mdl

s = m4.solve(tf=100.0, npoints=2000, outputs="x1 x2 x3".split())

def transformation(vars, t):
if t > 40.0:
return (vars[0] - 5.0, vars[1], vars[2])
else:
return (-5.0, vars[1], vars[2])

s.apply_transf(transformation)

s.plot()

title Rossler

X1' = X2 - X3
X2' = 0.36 * X2 - X1
X3' = X1 *X3 - 22.5 * X3 - 49.6 * X1 + 1117.8

init : (X1 = 19.0, X2 = 47, X3 = 50)
tf:200
~ x3 = X3 -50.0
~ x1 = X1 -18.0
~ x2 = X2 -50.0


## Example 5¶

Lorentz chaotic system: sensitivity to initial conditions

mdl = examples.models.lorentz.text
print mdl

ivs = {'init.x':(1.0, 1.01, 1.02)}
titles = ['$x(0)$ = %g' % v for v in ivs['init.x']]
s = m5.scan(ivs, tf=25.0, npoints=20000, outputs=['x'], titles=titles)

title Lorentz model: sensitivity to initial conditions
x' = 10*(y-x)
y' = x*(28-z)-y
z' = x*y - (8/3)*z
init: x = 1, y = 1, z = 1


## Example 6¶

CICR model again: parameter scanning

m = read_model("""
title Calcium Spikes
v0         = -> Ca, 1
v1         = -> Ca, k1*B*step(t, 1.0)
k1         = 7.3
B          = 0.4
export     = Ca ->  , 10 ..
leak       = CaComp -> Ca, 1 ..
!! Ca
v2         = Ca -> CaComp, 65 * Ca**2 / (1+Ca**2)
v3         = CaComp -> Ca, 500*CaComp**2/(CaComp**2+4) * Ca**4/(Ca**4 + 0.6561)
init       : (Ca = 0.1, CaComp = 0.63655)""")

import matplotlib as mpl
mpl.rcParams['figure.subplot.hspace']=0.4   #.2

Bstim = (0.0, 0.1, 0.2, 0.25, 0.28, 0.29, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.75, 0.8, 0.9, 1.0)
titles = ['$\\beta$ = %g' % B for B in Bstim]

s = m.scan({'B': Bstim}, tf=8.0, npoints=1000)
suptitlegend="Dynamics of cytosolic $Ca^{2+}$ as a function of stimulus"
s.plot(yrange=(0,1.5), legend=False, fig_size=(13.0,13.0), suptitlegend=suptitlegend)