# Format of ODE Files and Examples

## Contents

  Format for ODE files
  Miscellaneous functions
  Passive membrane
  Morris-Lecar equations
  Post-inhibitory rebound
  Hodgkin-Huxley
  Morris-Lecar with synapse
  Phase model
  Standard map
  Lorenz equations
  Unfolding of triple zero eigenvalue
  Cable equation
  Differential delay model
  A model with Markov processes
  Tyson's model with discontinuities
  Volterra integral equation

## ODE file format

Every ODE file consists of a series of lines that start with a keyword followed by numbers, names, and formulas or declare a named formula such as a differential equation or auxiliary quantity. Only the first letter of the keyword is important; e.g. the parser treats "parameter" and "punxatawney" exactly the same. The parser can understand lines up to 256 characters. You can use line continuation by adding a backslash character. The first line of the file cannot be a number (as this tells XPP that the file is in the old-style) but can be any other charcter or declaration. It is standard form to make the first line a comment which has the name of the file, but this is optional.

    # comments
     option  <filename>
    ...
    d<name>/dt=<formula>
    <name>'=<formula>
    ...
    <name>(t)=<formula>
    volt  <name>=<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
    ...
    @  <name>=<value>  ...
    done

- 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
  In this format the presence of the int operator tells XPP that the definition is for an integral equation and not a function definition. You can force XPP to treat the equation as a Volterra equation by strating the line with volt .
- Fixed or hidden values 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. They are just written as:

      zippy=f(x,t,...)

  Thus `zippy` could be used in seceral different places and yet is computed only once.
- Functions are easy to define. Just write like you would on paper:

      f(x,y) = x^2/(x^2+y^2)

  They can have up to 9 arguments.
- Initial data are written in two formats:

      x(0)=3
      y(0)=-2

  or

      init x=3,y=-2

  Initial data are optional, XPP sets them to zero by default and they can be changed within the program.
- Auxiliary quantities are expressions that depend on all of your dynamic variables which you want to keep track of. Energy is one such example. They are declared like fixed quantities, but are prefaced by `aux` . You have access to these within XPP and can plot them, etc.
- Markov variables are finite state variables whose dynamics are defined via the transition matrix that follows their declaration. They are declared by a statement that has the name of the process, the number of states and then the transition matrix. Initial values are as above. The entries of the transition 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.
- 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.
- Option files contain default values for such things as the step size for integration, the total amount of time to integrate, etc. Their use has been basically supplanted by the ability to declare such things within the ODE file using the @ operator.
- 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 .

# 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...
- shift(var,expr) This operator evaluates the expression expr converts it to an integer and then uses this to indirectly address a variable whose address is that of var plus the integer value of the expression. This is a way to imitate arrays in XPP. For example if you defined the sequence of 5 variables, u0,u1,u2,u3,u4 one right after another, then shift(u0,2) would return the value of u2.
- sum(ex1,ex2)of(ex3) is a way of summing up things. The expressions ex1>, are evaluated and their integer parts are used as the lower and upper limits of the sum. The index of the sum is i' so that you cannot have double sums since there is only one index. ex3 is the expression to be summed and will generally involve i' For example sum(1,10)of(i') will be evaluated to 55. Another example combines the sum with the shift operator. sum(0,4)of(shift(u0,i')) will sum up u0 and the next four variables that were defined after it.

### Internal Options

XPP has many many internal parameters that you can set from within the program and four parameters that can only be set before it is run. All of these internal parameters can be set from the ``Options'' files described above However, it is often useful to put the options right into the ODE file. NOTE: Any options defined in the ODE file overide all others such as the ones in the OPTIONS file. The format for changing the options is:

    @ name1=value1, name2=value2, ...

where name is one of the following and value is either an integer, floating point, or string. (All names can be upper or lower case). The first four options can only be set outside the program They are:

- MAXSTOR=`integer` sets the total number of time steps that will be kept in memory. The default is 5000. If you want to perform very long integrations change this to some large number.
- BACK= `{Black,White}` sets the background to black or white.
- SMALL=`fontname` where `fontname` is some font available to your X-server. This sets the ``small'' font which is used in the Data Browser and in some other windows.
- BIG=`fontname` sets the font for all the menus and popups.

The remaining options can be set from within the program. They are

- XP=name sets the name of the variable to plot on the x-axis. The default is `T` , the time-variable.
- YP=name sets the name of the variable on the y-axis.
- ZP=name sets the name of the variable on the z-axis (if the plot is 3D.)
- AXES=`{2,3}` determine whether a 2D or 3D plot will be displayed.
- TOTAL=value sets the total amount of time to integrate the equations (default is 20).
- DT=value sets the time step for the integrator (default is 0.05).
- NJMP=`integer` tells XPP how frequently to output the solution to the ODE. The default is 1, which means at each integration step.
- T0=value sets the starting time (default is 0).
- TRANS=value tells XPP to integrate until `T=TRANS` and then start plotting solutions (default is 0.)
- NMESH=`integer` sets the mesh size for computing nullclines (default is 40).
- METH=`{ discrete,euler,modeuler,rungekutta,adams,gear,volterra, backeul}` sets the integration method (see below; default is Runge-Kutta.)
- DTMIN=value sets the minimum allowable timestep for the Gear integrator.
- DTMAX=value sets the maximum allowable timestep for the Gear integrator
- TOLER=value sets the error tolerance for the Gear integrator.
- BOUND=value sets the maximum bound any plotted variable can reach in magnitude. If any plottable quantity exceeds this, the integrator will halt with a warning. The program will not stop however (default is 100.)
- DELAY=value sets the maximum delay allowed in the integration (default is 0.)
- XLO=value,YLO=value,XHI=value,YHI=value set the limits for two-dimensional plots (defaults are 0,-2,20,2 respectively.) Note that for three-dimensional plots, the plot is scaled to a cube with vertices that are +- 1 and this cube is rotated and projected onto the plane so setting these to +- 2 works well for 3D plots.
- XMAX=value, XMIN=value, YMAX=value, YMIN=value, ZMAX=value, ZMIN=value set the scaling for three-d plots.

# Examples

## The passive membrane

    # 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

The order of most declarations in the is unimportant since the program makes two passes through the file.

Return to tutorial

## Morris-Lecar Equations

    # 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 is:

    # 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

    # 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

    # 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:

    # 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 is:

    # 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 :

    # 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

    #map for triple-xero
    z(0)=2
    par a=2.93,b=-1.44,c=1.85
    z(t+1)= a+b*(z-c)^2
    done

The ODE file for computing the Liapunov exponent for the map is:

    # 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 is:

    # 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

Note that I have initialized V to be 1 which is the appropriate value to agree with the first boundary condition.

### 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 are:

    # 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.

Without further ado ...

    # 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.` Here it 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 convolution 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 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 ODE file 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