Tutorial - PyDSTool [Link]

[Link]

Site navigation:

ProjectOverview
CodeTopics
GettingStarted
Tutorial
UserDocumentation
ToolboxDocumentation
TechDocumentation

Elementary tutorial: A one-dimensional nonlinear

ODE (Calcium channel model)
See Tutorial_Calcium.py.

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), \quad \) (Ca)

where .

The system is specified using PyDSTool with

import PyDSTool as dst

import numpy as np
from matplotlib import pyplot as plt

# we must give a name

DSargs = [Link](name='Calcium channel model')
# parameters
[Link] = { 'vl': -60,
'vca': 120,
'i': 0,
'gl': 2,
'gca': 4,
'c': 20,
'v1': -1.2,
'v2': 18 }
# auxiliary helper function(s) -- function name: ([func signature], definition)
[Link] = {'minf': (['v'], '0.5 * (1 + tanh( (v-v1)/v2 ))') }
# rhs of the differential equation, including dummy variable w
[Link] = {'v': '( i + gl * (vl - v) - gca * minf(v) * (v-vca) )/c',
'w': 'v-w' }
# initial conditions
[Link] = {'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 cannot 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 Eq. (Ca) can be computed using a Generator instance:

[Link] = [0,30] # set the range of integration.

ode = [Link].Vode_ODEsystem(DSargs) # an instance of the 'Generator' class.
traj = [Link]('polarization') # integrate ODE
pts = [Link](dt=0.1) # Data for plotting

# PyPlot commands
[Link](pts['t'], pts['v'])
[Link]('time') # Axes labels
[Link]('voltage') # ...
[Link]([0,65]) # Range of the y axis
[Link]([Link]) # Figure title from model name
[Link]()

Depending on your local configuration of the Matplotlib interactive mode, the last command [Link]() might
not be necessary.

1 of 4 16/06/25, 8:18 pm
Tutorial - PyDSTool [Link]

The system described by Eq. (Ca) is bistable. This can be easily seen integrating trajectories with different
initial conditions:

[Link]() # Clear the figure

[Link](True) # Sequences of plot commands will not clear the existing figure
for i, v0 in enumerate([Link](-80,80,20)):
[Link]( ics = { 'v': v0 } ) # Initial condition
# Trajectories are called pol0, pol1, ...
# sample them on the fly to create Pointset tmp
tmp = [Link]('pol%3i' % i).sample() # or specify dt option to sample to sub-sample
[Link](tmp['t'], tmp['v'])
[Link]('time')
[Link]('voltage')
[Link]([Link] + ' multi ICs')
[Link]()

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

[Link](pars = {'i': -220} ) # Lower bound of the control parameter 'i'
[Link](ics = {'v': -170} ) # Close to one of the steady states present for i=-220

PC = [Link](ode) # Set up continuation class

PCargs = [Link](name='EQ1', type='EP-C') # 'EP-C' stands for Equilibrium Point Curve. The branch will be labeled 'EQ1'.
[Link] = ['i'] # control parameter(s) (it should be among those specified in [Link])
[Link] = 450 # The following 3 parameters are set after trial-and-error
[Link] = 2
[Link] = 1e-5
[Link] = 2e-2
[Link] = 'LP' # detect limit points / saddle-node bifurcations
[Link] = 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

2 of 4 16/06/25, 8:18 pm
Tutorial - PyDSTool [Link]

[Link](PCargs)
PC['EQ1'].forward()
[Link](['i','v'], stability=True, figure=3) # stable and unstable branches as solid and dashed curves, resp.

now consists of a "struct" data type that specifies the particular equilibrium curve we prepared the
PC['EQ1']
system for. The information of the equilibrium curve can be accessed via the info() method:
>>> PC['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 PC['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 = [Link](name='SN1', type='LP-C')

[Link] = 'EQ1:LP2'
[Link] = ['i', 'gca']
[Link] = 2
[Link] = ['CP']
[Link] = 200
[Link](PCargs)
PC['SN1'].forward()
PC['SN1'].backward()
PC['SN1'].display(['i','gca'], figure=4)

3 of 4 16/06/25, 8:18 pm
Tutorial - PyDSTool [Link]

4 of 4 16/06/25, 8:18 pm