## Tutorial

This tutorial illustrates how to use PyDSTool to analyze dynamical systems. It is inspired by the excellent XPP tutorial, from which it borrows some of the examples.

## A one-dimensional nonlinear ODE

A simple nonlinear model for the membrane voltage of a neuron is $C \frac{dV}{dt} = I + g_L (V_L - V) + g_{Ca} m(V)(V_{Ca} - V),$ where $m(V) = 0.5 ( 1 + \tanh[ (V-V_1) / V_2 ]).$

The system is specified using PyDSTool with

import PyDSTool
from pylab import plot, show, linspace, xlabel, ylabel

# we must give a name
DSargs = PyDSTool.args(name='Calcium')
# parameters
DSargs.pars = { 'vl': -60,
'vca': 120,
'i': 0,
'gl': 2,
'gca': 4,
'c': 20,
'v1': -1.2,
'v2': 18  }
# auxiliary helper function(s)
DSargs.fnspecs  = {'minf': (['v'], '0.5 * (1 + tanh( (v-v1)/v2 ))') }
# rhs of the differential equation, including dummy variable w
DSargs.varspecs = {'v': '( i + gl * (vl - v) - gca * minf(v) * (v-vca) )/c',
'w': 'v-w' }
# initial conditions
DSargs.ics      = {'v': 0, 'w': 0 }


The variable w is dummy, it merely tracks v. It is a necessary augmentation to permit the 2-parameter continuation of the limit points later on (which otherwise can't be done for a 1D dynamical system). This augmentation will be automatic and internal in a future version of PyCont.

### Integral curves

The solution of the dynamical system * can be computed using a Generator instance:

DSargs.tdomain = [0,40]                             # set the range of integration.
ode  = PyDSTool.Generator.Vode_ODEsystem(DSargs)    # an instance of the 'Generator' class.
traj = ode.compute('polarization')                  #
pd   = traj.sample()                                # Data for plotting
plot(pd['t'], pd['v'])
xlabel('time')                                      # Axes labels
ylabel('voltage')                                   # ...
ylim([0,65])                                        # Range of the y axis
show()


Depending on your local configuration, the last command 'show()' might not be necessary.

The system described by Equation * is bistable. This can be easily seen integrating trajectories with different initial conditions:

clf()                   # Clear the screen
hold(True)              # Sequences of plot commands will not clear the existing figures
for i, v0 in enumerate(linspace(-80,80,20)):
ode.set( ics = { 'v': v0 } )                       # Initial condition
tmp = ode.compute('pol%3f' % i).sample()           # Trajectories are called pol0, pol1, ...
plot(tmp['t'], tmp['v'])
xlabel('time')
ylabel('voltage')
show()


### Bifurcation diagrams

To see how the fixed points depend on the parameters we will plot bifurcation diagrams using the continuation class (ContClass). We start with the diagram that shows the equilibrium voltage v as a function of the input i.

# Prepare the system to start close to a steady state
ode.set(pars = {'i': -220} )       # Lower bound of the control parameter 'i'
ode.set(ics =  {'v': -170} )       # Close to one of the steady states present for i=-220

PyCont = PyDSTool.ContClass(ode)                 # Set up continuation class

PCargs = PyDSTool.args(name='EQ1', type='EP-C')  # 'EP-C' stands for Equilibrium Point Curve. The branch will be labeled 'EQ1'.
PCargs.freepars     = ['i']                      # control parameter(s) (it should be among those specified in DSargs.pars)
PCargs.MaxNumPoints = 450                        # The following 3 parameters are set after trial-and-error
PCargs.MaxStepSize  = 2
PCargs.MinStepSize  = 1e-5
PCargs.StepSize     = 2e-2
PCargs.LocBifPoints = 'LP'                       # detect limit points / saddle-node bifurcations
PCargs.SaveEigen    = True                       # to tell unstable from stable branches


The LocBifPoints attribute tells PyCont what type of bifurcation should be tracked (see PyCont for details). In this example we specify that only saddle-node bifurcations ('LP') should be detected. The SaveEigen attribute is a boolean variable that determines whether or not the eigenvalues of the equilibrium points should be saved along the curve. We set this attribute to True because we want to know the stability along equilibrium curve. Once the continuation class is set up, we can compute the bifurcation diagram

PyCont.newCurve(PCargs)
PyCont['EQ1'].forward()
clf()
PyCont.display(['i','v'], stability=True)        # stable and unstable branches as solid and dashed curves, resp.


PyCont['EQ1'] now consists of a "struct" data type that specifies the particular equilibrium curve we prepared the system for. The information of the equilibrium curve can be accessed via the info() method:

>>> PyCont['EQ1'].info()
PyCont curve EQ1 (type EP-C)
Using model: Calcium_model

Model Info
----------
Variables : v
Parameters: gca, vca, c, i, vl, v1, v2, gl

Continuation Parameters
-----------------------
name          = EQ1
freepars      = ['i']
auxpars       = []
MaxNumPoints  = 450
MaxCorrIters  = 5
MaxTestIters  = 10
MaxStepSize   = 2
MinStepSize   = 1e-05
StepSize      = 2
VarTol        = 1e-06
FuncTol       = 1e-06
TestTol       = 0.0001
LocBifPoints  = ['LP']
...

Special Points
--------------
P1, P2, LP1, LP2



We can obtain detailed information about a particular special point calling the getSpecialPoint method. For instance, limit point LP2 has the following properties:

>>> print PyCont['EQ1'].getSpecialPoint('LP2')
i:  -210.477955042
v:  15.4485534278
..


We now want to know the location of the limit points change as we vary the calcium conductance, i.e., the parameter gca. We start from one of the limit points, say LP2,

PCargs = PyDSTool.args(name='SN1', type='LP-C')
PCargs.initpoint    = 'EQ1:LP2'
PCargs.freepars     = ['i', 'gca']
PCargs.MaxStepSize  = 2
PCargs.LocBifPoints = ['CP']
PCargs.MaxNumPoints = 200
PyCont.newCurve(PCargs)
PyCont['SN1'].forward()
PyCont['SN1'].backward()
PyCont['SN1'].display(['i','gca'])


PyDSTool source code is hosted by:

last edited 2010-03-26 22:47:37 by RobClewley