Format of ODE Files and Examples

Every ODE file consists of a series of lines that start with a keyword followed by numbers, names, and formulas. Only the first letter of the keyword is important; e.g. the parser treats "variable" and "vavoom" exactly the same. The parser can understand lines up to 256 characters. You can not use line continuation in the old style parser. The new style parser is much more intuitive and flexible. You should read this but use that.


<number>
change <filename>
variable <name>=<value>, ... 
...
fixed <name>,...
...
aux <name>,...
...
markov <name> <value> <nstates>
...
wiener <name>,...
...
global sign {condition} {name1=form1;name2=form2;...}
...
parameter <name>=<value>, ...
...
number <name>=<value>
...
table <name> <filename>
...
table <name> % <npts> <tlo> <thi> <function(t)>
...
user <name> <nargs> <expression>
...
kernel <name> <mu> <expression>
...
ode <expression>
...
integral <expression>
...
ode <expression>   (for fixed and auxilliary quantities)
...
bdry <expression>
...
r <name>
{t01} {t02} ... {t0k-1}
{t10} ...
...
{tk-1,0} ...

# Comments      

done

Every ODE file must have a number as its first line. This number is the total number of dependent dynamic quantities (such as variables and Markov variables) and auxiliary quantities. If the first line is not a number the new style for ODE files is assumed.
change This tells XPP to use the file specified instead of "default.opt." These files contain default values for such things as the step size for integration, the total amount of time to integrate, etc.
variable This declares the name of a variable which satisfies some differential, integral, or discrete equation. Initial values are optional.
fix Fixed variables are quantities that you define but do not satisfy differential equations and which are not accessible to the user. They are generally expressions that are used in several places in the dynamical system so that you avoid computing the same quantity several times. You just list their names at the beginning of the file.
markov Markov variables are finite state variables whose dynamics are defined via the transition matrix. A single line is required for each one. The name, the initial state and the number of states are separated by spaces. They must be declared after all fixed and regular variables are defined.
auxiliary Auxiliary quantities are expressions that depend on all of your dynamic variables which you want to keep track of. Energy is one such example. Their names are just listed separated by spaces or commas.
wiener Wiener parameters are normally distributed numbers with zero mean and unit standard deviation. They are useful in stochastic simulations since they automatically scale with change in the integration time step. Their names are listed separated by commas or spaces.
global Global flags are expressions that signal events when they change sign, from less than to greater than zero if sign=1 , greater than to less than if sign=-1 or eithet way if sign=0. The condition should be delimited by braces {} The events are of the form variable=expression , are delimited by braces, and separated by semicolons. When the condition occurs all the variables in the event set are changed possibly discontinuously. This lets you simulate "Delta" functions.
parameter This is how you declare parameters that you may want to change. The values are optional; they are otherwise set to zero.
number Number declarations are like parameter declarations, except that they cannot be changed within the program and do not appear in the parameter window.
table The Table declaration allows the user to specify a function of 1 variable in terms of a lookup table which uses linear interpolation. The name of the function follows the declaration and this is followed by (i) a filename (ii) or a function of "t". The file is an ascii readable file in this format:
```
number of values
xlo
xhi
y(xlo)
.
.
.
y(xhi)
```
The function form of the table require the percent symbol followed by the number of points, tlo, thi and the function of t.
User functions XPP allows you to define your own functions of up to 9 variables. The variables are named arg1,arg2,...,arg9. After the User declaration, follow with the name of the function, the number of arguments, and the definition.
kernel Declare a kernal for a Volterra integral equation
ode The ode declaration is how you assign the right-hand sides to differential equations, fixed variables, and auxiliary variables. The order is important . The differential equations must be defined first in the order of the variable you have declared. Then the fixed quantities. Finish with the auxiliary variables.
integral The Integral declaration is like the Ode declaration but tells the parser that the equation is not a differential equation, but instead, an integral equation.
bndry The boundary declaration tells XPP about boundary values. Since a boundary value problem is defined on an interval [a,b] XPP uses the variable names to mean the value at a and primed variable names to mean the value at b. The program tries to set the expression to zero. Thus x implies that x(a)=0 and x'-2 implies that x(b)-2=0 .
R This declaration defines the transition matrix for the named markov variable. The entries of the matrix are delimited by braces. The first row represents the transition probability of switching from state 1 to some other state. The entries of the matrix can be any expression. The diagonal entries are ignored.

Miscellaneous functions

heav(x) the step function.
sign(x) is 1 for x>0, -1 for x<0, and 0 for x=0.
ran(arg) produces a uniformly distributed random number between 0 and arg.
besselj, bessely take two arguments, n,x and return respectively, J(n,x) and Y(n,x).
erf(x),erfc(x) are the Error function and the complementary error function.
normal(arg1,arg2) produces a normally distributed random number with mean arg1 and variance arg2.
max(arg1,arg2) produces the maximum of the two arguments and min(arg1,arg2) , the minimum of them.
if(expr1)then(expr2)else(expr3) evaluates expr1 If it is nonzero it evaluates to expr2 otherwise it is expr3 .
delay(var,expr) returns variable var delayed by the result of evaluating . In order to use the delay you must inform the program of the maximal possible delay so it can allocate storage.
ceil(arg),flr(arg) are the integer parts of arg returning the smallest integer greater than and the largest integer less than arg .
t is the current time in the integration of the differential equation.
pi is 3.1415926...

New style parser

The newstyle parser lets you use the standard UNIX line continuation \ so that you do not have to make the lines very long.

The newstyle format for the parser has the following format:

 option  <filename>
...
d<name>/dt=<formula>
<name>'=<formula>
...
<name>(t)=<formula>
...
<name>(t+1)=<formula>
...
<name>=<formula>
...
<name>(<x1>,<x2>,...,<xn>)=<formula>
...
init <name>=<value>,...
...
<name>(0)=<value>
...
aux <name>=<formula>
...
bdry <expression>
...
wiener <name>,...
...
parameter <name>=<value>, ...
...
number <name>=<value>
...
table <name> <filename>
...
table <name> % <npts> <xlo> <xhi> <function(t)>
...
global sign {condition} {name1=form1;...}
...
markov <name> <nstates>
{t01} {t02} ... {t0k-1}
{t10} ...
...
{tk-1,0} ...
...
# comments
...
done

The big difference is that the name of the variable, auxiliary quantity, fixed quantity, and Markov variable are kept with their right-hand sides. Thus there is no o,i,r stuff necessary anymore and the fixed,variable,kernel declarations are gone. Auxiliary variables always have user defined names and functions are put in as you would write them. (No longer is it necessary to use arg1,arg2, etc as the names of the arguments.) For variables that will satisfy differential equations, maps, or integral equations, you write the equation in the most obvious fashion.

Ordinary differential equations are written as
```
dx/dt=F(x,y,...)
```
or
```
x'=F(x,y,...)
```
Similarly, difference equations can be written as:
```
x(t+1)=F(x,y,...)
```
Volterra integral equations are written in one of the following 4 formats:
```
x(t) =  ...int{K(u,t,t')}...
x(t) =  ...int[q]{K(u,t,t')}...
x(t) =  ...int{k(t)#u}...
x(t) =  ...int[q]{k(t)#u}...
```
with the following meanings
- General integral
- Integral with singularity of the form 1/t^q
- Convolution with k(t)
- Convolution with k(t)/t^q
Fixed values are just written as:
```
zippy=f(x,t,...)
```
Functions are much easier to define. Just write is like you would on paper:
```
f(x,y) = x^2/(x^2+y^2)
```
Initial data can be wriiten in two formats:
```
x(0)=3
y(0)=-2
```
or
```
init x=3,y=-2
```
Auxiliary quantities are wriiten like fixed quantities, but are prefaced by aux.
Markov processes are defined by a statement that has the name of the process, the number of states and then the transition matrix. Initial values are as above.
Option files are now declared with the option statement instead of the change statement.

See the above section for details on the rest of the declarations.

Examples

The passive membrane

Old Style:

New Style:

For completeness, I also have included the new style for the ODE files following the old style. So for the above the format is:

# passive membrane model

# The parameters:
parameter c=20,gl=2,vl=-60,i=2

# the differential equation:
dv/dt=(gl*(vl-v)+i)/c

# initial data:
v(0)=0
# The initial data is not needed here as the default is 0

# tell it you are done  
done

It is pretty much self-explanatory. Initial data are declared like:

name(0) = value

where "name" is the name of a variable which satisfies a differential equation. This is fine for one or two variable systems, but if you have many, then it is a pain to write each one on a line. Instead, you can write:

init name1=value1, name2=value2, ...

Differential equations in the new format are declared as

name' = right-hand-side
# OR
dname/dt = right-hand-side

You also do not need to tell the program how many equations. Finally, the order of most declarations in the new style is unimportant since the program makes two passes through the file.

Return to tutorial

Morris-Lecar Equations

Old Style:

I have added some comments

4
# Morris-Lecar Equations with ICa and IK tracked
# name and initialize variables:
v v=-60.899,w=.014873
aux ICa,IK
# name and initialize parameters:
param vk=-84,vl=-60,vca=120
param i=0,gk=8,gl=2, gca=4, c=20
# type II -- v3=2,v4=30,phi=.04 -- like HH 
# type I -- v3=12,v4=17.4,phi=.06666667
param v1=-1.2,v2=18,v3=2,v4=30,phi=.04
# define the functions
user minf 1 .5*(1+tanh((arg1-v1)/v2))
user winf 1 .5*(1+tanh((arg1-v3)/v4))
user lamw 1 phi*cosh((arg1-v3)/(2*v4))
# The right-hand sides
odev  (i+gl*(vl-v)+gk*w*(vk-v)+gca*minf(v)*(vca-v))/c
odew lamw(v)*(winf(v)-w)
# the auxiliary variables
oICa gca*minf(V)*(V-Vca)
oIK  gk*w*(V-VK)
done

New Style:

# Morris-Lecar reduced model 
dv/dt=(i+gl*(vl-v)+gk*w*(vk-v)+gca*minf(v)*(vca-v))/c
dw/dt=lamw(v)*(winf(v)-w)
# where
minf(v)=.5*(1+tanh((v-v1)/v2))
winf(v)=.5*(1+tanh((v-v3)/v4))
lamw(v)=phi*cosh((v-v3)/(2*v4))
#
param vk=-84,vl=-60,vca=120
param i=0,gk=8,gl=2, gca=4, c=20
param v1=-1.2,v2=18,v3=2,v4=30,phi=.04
# for type II dynamics, use v3=2,v4=30,phi=.04
# for type I dynamics, use v3=12,v4=17,phi=.06666667
v(0)=-60.899
w(0)=0.014873
# track some currents
aux Ica=gca*minf(V)*(V-Vca)
aux Ik=gk*w*(V-Vk)
done

Return to tutorial

Post-inhibitory Rebound

This model is the simplest neural oscillator. There is only one active channel. The equations are due to Wang and Rinzel. The equations include the activation of the fast GABA_A and slow GABA_B synapses.

A good set of parameters is: gca=1,vk=-90,fh=2,vca=120, gl=.1,vl=-70,i=0, gsyna=0,gsynb=0,vsyna=-85,vsynb=-90,tvsyn=-40, fka=2,rka=.08. Try looking at the phase-plane with -75 < V < 20 and -.125 < h < .5. By letting the current be slightly negative, try and get this baby to oscillate. Look at the nullclines as a function of the current.

Hint Integrate it for about 200 msec with DeltaT=.1.

If you still don't think you can create the ODE file yourself, click here!

Go to part 4

The Hodgkin-Huxley Equations

The Hodgkin-Huxley equations are the classic model for action potential propagation in the squid giant axon. They represent a patch of membrane with three channels, sodium, potassium, and a leak. The equations are:

You should explore them as the current, I increases. Try to find the Hopf bifurcation. Show that there is a region where there is bistability between a periodic and fixed point. This bistability was the basis for an experiment done by Rita Gutman in the 70's.

The ODE file in the old style is:

4
v v,m,h,n
p vna=50,vk=-77,vl=-54.4,gna=120,gk=36,gl=.3,c=1,i=0
u am 1 .1*(arg1+40)/(1-exp(-(arg1+40)/10))
u bm 1 4*exp(-(arg1+65)/18)
u ah 1 .07*exp(-(arg1+65)/20)
u bh 1 1/(1+exp(-(arg1+35)/10))
u an 1 .01*(arg1+55)/(1-exp(-(arg1+55)/10))
u bn 1 .125*exp(-(arg1+65)/80)
oV (I - gna*h*(v-vna)*m^3-gk*(v-vk)*n^4-gl*(v-vl))/c
om am(v)*(1-m)-bm(v)*m
oh ah(v)*(1-h)-bh(v)*h
on an(v)*(1-n)-bn(v)*n
d

and in the new style:

# Hodgkin-huxley equations
dv/dt=(I - gna*h*(v-vna)*m^3-gk*(v-vk)*n^4-gl*(v-vl))/c
dm/dt= am(v)*(1-m)-bm(v)*m
dh/dt=ah(v)*(1-h)-bh(v)*h
dn/dt=an(v)*(1-n)-bn(v)*n
par vna=50,vk=-77,vl=-54.4,gna=120,gk=36,gl=.3,c=1,i=0
am(v) =  .1*(v+40)/(1-exp(-(v+40)/10))
bm(v) =  4*exp(-(v+65)/18)
ah(v) =  .07*exp(-(v+65)/20)
bh(v) =  1/(1+exp(-(v+35)/10))
an(v) =  .01*(v+55)/(1-exp(-(v+55)/10))
bn(v) =  .125*exp(-(v+65)/80)
init v=-65,m=.052,h=.596,n=.317
done

Back to the tutorial

Morris-Lecar with synapse

In the old style

3
v v=-24.22,w=.305,s=.056
param vk=-84,vl=-60,vca=120
param i=80,gk=8,gl=2, gca=4, c=20
# type II -- v3=2,v4=30,phi=.04 -- like HH 
# type I -- v3=12,v4=17.4,phi=.06666667
param v1=-1.2,v2=18,v3=12,v4=17.4,phi=.06666667
param vsyn=120,vt=20,vs=2,alpha=1,beta=.25
user minf 1 .5*(1+tanh((arg1-v1)/v2))
user winf 1 .5*(1+tanh((arg1-v3)/v4))
user lamw 1 phi*cosh((arg1-v3)/(2*v4))
user k 1 1/(1+exp(-(arg1-vt)/vs))
odev  (i+gl*(vl-v)+gk*w*(vk-v)+gca*minf(v)*(vca-v))/c
odew lamw(v)*(winf(v)-w)
# synapse
odes alpha*k(v)*(1-s)-beta*s
done

and in the new style:

# Morris-Lecar reduced model 
dv/dt=(i+gl*(vl-v)+gk*w*(vk-v)+gca*minf(v)*(vca-v))/c
dw/dt=lamw(v)*(winf(v)-w)
ds/dt=alpha*k(v)*(1-s)-beta*s
minf(v)=.5*(1+tanh((v-v1)/v2))
winf(v)=.5*(1+tanh((v-v3)/v4))
lamw(v)=phi*cosh((v-v3)/(2*v4))
k(v)=1/(1+exp(-(v-vt)/vs))
param vk=-84,vl=-60,vca=120
param i=80,gk=8,gl=2, gca=4, c=20
param v1=-1.2,v2=18,v3=12,v4=17.4,phi=.066666667
param vsyn=120,vt=20,vs=2,alpha=1,beta=.25
# for type II dynamics, use v3=2,v4=30,phi=.04
# for type I dynamics, use v3=12,v4=17.4,phi=.06666667
v(0)=-24.22
w(0)=.305
s(0)=.056
done

Back to the salt mines ...

A simple phase-difference model

In the old style

2
v x1,x2
p a0=.667,delta=0,a1=-.98,b1=.8174,g=.25
u h 1 a0+a1*cos(arg1)+b1*sin(arg1)
o delta+g*h(x2-x1)
o g*h(x1-x2)
d

and in the new style

# phase model based on ML example
p a0=.667,delta=0,a1=-.98,b1=.8174,g=.25
h(phi)=a0+a1*cos(phi)+b1*sin(phi)
x1'=delta+g*h(x2-x1)
x2'=g*h(x1-x2)
d

Return to sender

Standard map

This discrete dynamical system has the oldstyle format:

1
v y
par g=.35,delta=.5
o y+delta-g*sin(y)
d

and in the newstyle:

# standard map
y(n+1)=y+delta-g*sin(y)
par g=.35,delta=.5
done

Lorenz equations

The Lorenz equations are a three dimensional dynamical system representing turbulence in the atmosphere:

where r=27,s=10,b=8/3. The ODE file in the old style format is:

3
# the famous Lorenz equation
v x=1,y=2,z
p r=27,s=10,b=2.66666
o s*(-x+y)
o r*x-y-x*z
o -b*z+x*y
d

and in the new style

# the famous Lorenz equation
init x=1,y=2
p r=27,s=10,b=2.66666
x'= s*(-x+y)
y'= r*x-y-x*z
z'=  -b*z+x*y
d

Unfolding of the triple zero

The ODE file in the old style is:

3
v x,y,z
p mu=1,nu=2,gamma=1,q=2
o y-mu*x
o z-nu*x
o q*x*x-gamma*x
d

and in the new style

# triple zero unfolding
x'=y-mu*x
y'=z-nu*x
z'=q*x*x-gamma*x
p mu=1,nu=2,gamma=1,q=2
done

The ODE file for the approximate map for this system is

1
v z=2
p a=2.93,b=-1.44,c=1.85
o a+b*(z-c)^2
d

The ODE file for computing the Liapunov exponent for the map in the old style is:

3
v z=2,zp=0
aux lexp
p a=2.93,b=-1.44,c=1.85
o a+b*(z-c)^2
o zp+log(abs(2*b*(z-c)))
o zp/(t+1)
d

and in the new style

# Liapunov exponent
init z=2,zp=0
aux lexp=zp/(t=1)
par a=2.93,b=-1.44,c=1.85
z(n+1)= a+b*(z-c)^2
zp(n+1)= zp+log(abs(2*b*(z-c)))
d

Other examples with some commentary

A boundary value problem: steady state cable equation with leak

In studying a dendrite, one is often interested in the steady-state voltage distribution. Consider a case where the dendrite is held at V=V_0 at x=0 and has a leaky boundary condition at x=1. Then the steady state cable equations are:

Hint Solve this equation by setting the total integration time to 1 and then using the boundary-value solver command (B) (S). Try ranging from b=0 to b=10 and plotting the results in different colors.

Since this is a second order equation and XPP can only handle first order, we write it as a system of two first order equations in the file which in the old style is:

2
# define the variable, v and its derivative vx
v v=1,vx=0
# define the 4 parameters, a,b,v0,lambda^2
p a=1,b=0,v0=1,lam2=1
# now do the right-hand sides
o vx
o v/lam2
#
# and finally, boundary conditions
# First we want V(0)-V0=0
b v-v0
#
# We also want aV(1)+bVX(1)=0
b a*v'+b*vx'
# Note that the primes tell XPP to evaluate at the right endpoint
done

Note that I have initialized V to be 1 which is the appropriate value to agree with the first boundary condition. The equations in the new style are:

# Linear Cable Model

# define the 4 parameters, a,b,v0,lambda^2
p a=1,b=0,v0=1,lam2=1
# now do the right-hand sides
v'=vx
vx'=v/lam2

# The initial data
v(0)=1
vx(0)=0

# and finally, boundary conditions
# First we want V(0)-V0=0
b v-v0
#
# We also want aV(1)+bVX(1)=0
b a*v'+b*vx'
# Note that the primes tell XPP to evaluate at the right endpoint
d

Delayed inhibitory feedback net

A simple way to get oscillatory behavior is to introduce delayed inhibition into a neural network. The equation is

Hint Set the maximal delay to 10 by clicking (U) (E), setting it, and returning to the main menu. Then integrate the equations varying tau which is the delay. For what value of delay is the rest state unstable?

The equations in the old style are:

1
# declare the scalar x
v x=1
# declare all the parameters, initializing the delay to 3
p tau=3,b=4.8,a=4,p=-.8
# define the nonlinearity
u f 1 1/(1+exp(-arg1))
# define the right-hand sides; delaying x by tau
o -x + f(a*x-b*delay(x,tau)+p)
# done
d

and in the new style :

# delayed feedback

# declare all the parameters, initializing the delay to 3
p tau=3,b=4.8,a=4,p=-.8
# define the nonlinearity
f(x)=1/(1+exp(-x))
# define the right-hand sides; delaying x by tau
dx/dt = -x + f(a*x-b*delay(x,tau)+p)
x(0)=1
# done
d

A population problem with random mutation

This example illustrates the use of Markov variables and is due to Tom Kepler. There are two variables, x1,x2 and a Markov state variable, z. The equations are:

where the Markov variable z has two states and a transition matrix dependent on x1 given by:

Hint Set the total integration time to 40, window the X-axis between 0 and 40, and add a graph of x2 along with x1. Integrate this several times to see that x2 spontaneously arises once x1 is large enough and then like all really cool mutants takes over and kills x1.

The old style format is:

3
v x1=1.e-4,x2=1.e-4
m z 0 2
p eps=.1,a=1
o x1*(1-x1-x2)
o z*(a*x2*(1-x1-x2)+eps*x1)
r z
{0} {eps*x1}
{0} {0}
d

and the new style format is:

# Kepler model

# another way to do initial data
init x1=1.e-4,x2=1.e-4
p eps=.1,a=1
x1' = x1*(1-x1-x2)
x2'= z*(a*x2*(1-x1-x2)+eps*x1)
# the markov variable and its transition matrix
markov z 2
{0} {eps*x1}
{0} {0}
d

A system with discontinuities

John Tyson describes a model which involves cell growth:

This is a normal looking model with exponential growth of the mass. However, if u decreases through 0.2, then the ``cell'' divides in half. That is the mass is set to half of its value.

Hint Set the total integration time to 800 and nOut to 10, integrate the equations, and then plot the mass versus t .

We want to flag the condition u=.2 with u decreasing through .2 so we want the sign=-1. The source code in the new style (since older versions do not support global there is no sense using the old style, although it could be done) is:

# tyson
i u=.0075,v=.48,m=1
p k1=.015,k4=200,k6=2,a=.0001,b=.005
u'=  k4*(v-u)*(a+u^2) - k6*u
v'= k1*m - k6*u
m'= b*m
global -1 {u-.2} {m=.5*m} 
done

A simple Volterra integral equation

Suppose we want to solve:

Note that it is a convoltuion equation, so XPP can take advantage of that and speed up the problem considerably. I show you the ODE file with and without the convoltuion speed up. The closed form solution to this equation is u(t)=sin(t).

Hint Plot the solution along with the auxiliary variable, utrue to see that the solutions are the same.

The source without using the convolution speed up in the old style is:

2
# Volterra equation with actual solution
v u
a utrue
i sin(t)+.5*cos(t)-.5*t*exp(-t)-.5*exp(-t)+int{(t-t')*exp(-(t-t'))*u}
o sin(t)
d

and in the new style is:

# Volterra equation with actual solution
u(t) = sin(t)+.5*cos(t)-.5*t*exp(-t)-.5*exp(-t)+int{(t-t')*exp(-(t-t'))*u}
a utrue=sin(t)
d

A much faster version takes advantage of the convolution structure. The old style ODE file is:

2
# Faster version of testvol.ode
v u
a utrue
i sin(t)+.5*cos(t)-.5*t*exp(-t)-.5*exp(-t)+int{t*exp(-t)#u}
o sin(t)
d

and in the new style is:

# Faster version of testvol.ode
u(t)=sin(t)+.5*cos(t)-.5*t*exp(-t)-.5*exp(-t)+int{t*exp(-t)#u}
a utrue=sin(t)
done

Have a nice day!

References

Table of Contents
Main Menu Items
ODE Files and Examples
Numerics Menu
File Menu
Freeze Menu
AUTO Menu
Data Browser
I/O and Hardcopy
XPP Basics
Nonlinear ODEs
Two-dimensions
Three-dimensions and Beyond
Phase Equations
Chaos