Robert A McDougal: XPP Tutorial (2024)

Getting Started with XPPAUT

XPPAUT is aprogram by Bard Ermentrout of the University ofPittsburgh for solving a variety of dynamical systems problems on thecomputer. His official tutorial is located here, butwe will only concern ourselves with the basics.

Outline of Contents

  • Plotting Solutions to Ordinary Differential Equations
    • A Simple, First Order ODE
    • Changing the Window
    • Setting the Initial Conditions
    • Saving Images as Postscript Files
    • Displaying Numerical Values of Solution
    • A Second Order ODE
    • Converting a High Order ODE to a First Order System
    • Change Variables Plotted
    • Add a Curve to the Graph
    • PIR in Mutually Coupled FitzHugh-Nagumo Neurons
    • Parameters
    • Functions
    • Logical Expressions and Split Functions
    • Arrays
    • Configuring XPPAUT from the .ODE File
  • Delay Equations
    • Auxillary Functions
  • Phase Planes
    • Use the Mouse to Select Initial Conditions
    • Flow
    • Nullclines
    • Fixed Points
  • Bifurcation Diagrams

Plotting Solutions to Ordinary DifferentialEquations

A Simple, First Order ODE

We start by concerning ourselves with one of the simplest Initial ValueProblems that we can think of,

x' = 1
x(0) = −1

This problem clearly has solution x(t)=−1+t.We graph it with XPPAUT.First, let's create a text file sample1.ode:

# sample1.ode
# Solves the initial value problem dx/dt=1, x(0)=-1
dx/dt=1
x(0)=-1
done

The first two lines are comments and are ignored by XPP. In the thirdline, we specify the differential equation. The fourth line gives the initialvalue condition, and the fifth line tells XPP that we are done stating thesystem.

Warning: Do not be deceived by thestatement x(0)=-1. You cannot say x(1)=-1 if you instead wanted to start integrating fromt=0. The general way to state initial conditions is to use the form init x=-1 or just i x=-1.

We have the computer solve our ODE by typing xppautsample1.ode at the command line. A window appears:

Click Initialconds and then (O)ld to plot using the initial conditions wespecified in the .ode file. (Note: In XPPAUT, youcan also just type the capitalized letter in a menu choice instead ofclicking it. In our case then, we would just type io to graph.)

We have a graph, but the axes need adjusting. Click Viewaxes, then 2Dto display the axes menu.

Here we can specify which variables we want to plot, the range of x and yvalues, and labels for the axes. In the image above, I have changedthe y values to run from −5 to 5 and the x values to run from 0 to 5.I also put labels on the axes.

To change the initial condition, click on the cyan colored button ICs in the upper left-hand corner. A new windowappears.

Change your initial condition and click Go. A new line is drawn on your graph. To clear thegraph before drawing a new trajectory, click Erase (or type e)before clicking Go. I for example chose toadd a trajectory starting with x = −3.

Warning: It is important to note thatXPPAUT only stores your graph data if youexplicitly tell it to by selecting Graphicstuff and then (F)reeze, which storesthe last trajectory plotted. If you do not make that selection, the oldtrajectory will be erased if the window is redrawn, it will not be saved ifyou try to save the image, etc...

To save our trajectories as a postscript file, select Graphic stuff and then (P)ostscript. For now, hit Ok on the PostscriptParameters window, enter a filename and then click Ok to save.

Plotting trajectories can help you get an intuitive sense of what is goingon, but in many situations you want to know specific function values, eitherfor their own sake or so you can do further analysis on the problem inanother program. To do so, click on the cyan colored button Data and a newwindow will appear.

Note that the time values are spaced every . 05. This is the default. To change this behavior, select nUmerics and then Dt. Options for scrolling and writing the data todisc are available at the top of the window.

A Second Order ODE

Consider the ODE

x'' = −x
x(0) = 0
x'(0) = 1

which has solution x(t) = sin(t). XPPAUT only works with first order equations orsystems of equations, not second order. Fortunately by making the change ofvariables

x1 = x
x2 = x'

we can turn any second order ODE into a system of two first order ODEs.Similarly, any higher order ODE can also be reduced to a system of firstorder ODEs. In our case, for example, the system becomes:

x1' = x2
x2' = −x1
x1(0) = 0
x2(0) = 1

Our system then has solution x(t)=x1(t)=sin(t)and x2(t)=cos(t).We enter the problem into XPPAUT assample2.ode:

# sample2.ode
# solves the system form of x''=-x, x(0)=0, x'(0)=1
dx1/dt=x2
dx2/dt=-x1
x1(0)=0
x2(0)=1
done

As before, we run it by typing xppautsample2.ode at the command prompt, clicking on the cyan button ICs and selecting Go. We then have the following window.

Note that it plotted X1 vs t. To see X2 vs t instead, select Viewaxes then 2Dand change Y-axis to X2. (We could even plot X1 vs X2 -- which in our casewould give us a circle. See the Phase Planesection for more.)

Click Ok and the sine wave should bereplaced with a cosine wave.

We may prefer to have both X1 and X2 plotted at the same time. For example, in a laterexample we consider the interactions of two mutually coupled neurons, and wewould like to see both of their electrical outputs simultaneously.Fortunately, XPPAUT gives us an easy way toadd a curve to the current graph. Select Graphicstuff, then (A)dd curve; a new dialogbox appears.

Here I have changed Y-axis to x1 (case does not matter in XPPAUT) and Colorto 1. (The Z-axisoption does not matter since we are only working with a two-dimensionalimage.) Click Ok and the curve will be added.(The new curve will be red since we set Colorto 1; other color choices are available.)

PIR in Mutually Coupled FitzHugh-NagumoNeurons

We now attempt sample3.ode, a morecomplicated example that introduces a number of new concepts.

# sample3.ode
# Post Inhibitory rebound in mutually coupled
# FitzHugh-Nagumo cells.

# gsyn=2 gives suppressed, gsyn=.5 gives pir

# Parameter values. It's always better to give constants
# as parameters so you can change them without exiting
# XPP.
p vSpike=-1, eps=.2, gsyn=.5, theta=-.2,scale=10
p wshift=3.8, vshift=0, scale2=10
p alpha=1, beta=.1, vstart=0, wstart=0, vsyn=-2

# Usually we talk of the FitzHugh-Nagumo model with
# f and g somewhat arbitrary, but to do numerics we
# need to pick specific functions.
f(a,b)=-b-15*(a-1)*a*(a+1)
g(a,b)=-b+scale*tanh(a*scale2)+wshift

# the slow system (equation) and the synapse obey
# the same rules for both neurons, so we use this array
# notation. The first line denotes the range of j (always
# j) and v[j] denotes the j'th v value, etc...
%[1..2]
w[j]'=eps*g(v[j],w[j])
s[j]'=alpha*(1-s[j])*(v[j]>theta)-beta*s[j]
%

# cell 1 receives input from 2 and 2 from 1, so we write
# these seperately since they are different
v1'=f(v1,w1)-gsyn*(v1-vsyn)*s2
v2'=f(v2,w2)-gsyn*(v2-vsyn)*s1

# initial conditions. Since I don't specify s values,
# they start at 0.
i v1=0, w1=-5, v2=-1, w2=-2

# configure XPP's options.
@ total=100, xlo=0, xhi=100, ylo=-1.5, yhi=1.5
@ nplot=2, xp[1..2]=t, yp[j]=v[j]
done

Integrating using the default settings gives us the following graph. Notehow the membrane potential of each neuron drops a little when the other isfiring; the drop is because the neurons inhibit each other.

Parameters

Sometimes our system is critically dependent on the particular numbers inour equations. (Bifurcation diagrams helpus understand this dependency.) When you first attempt to model a problem,you may not know the constants. Using parameters (the lines beginning withthe letter p) we can change the constants withinXPPAUT, by simply clicking on the cyancolored button labelled Param.

Exercise. Vary the parameters to seewhat other behaviors you can get out of the system. Can you make the neuronsfire together? Can you make it so that only one neuron is firing?

Functions

Notice the lines

f(a,b)=-b-15*(a-1)*a*(a+1)
g(a,b)=-b+scale*tanh(a*scale2)+wshift

They define two functions, f and g, that I can use elsewhere in my code. The ability todefine functions helps with readablility and maintainability of your .ODEfiles.

Logical Expressions and SplitFunctions

As a first approximation, a synapse only experiences a buildup ofneurotransmitters when the membrane potential of the presynaptic cell isabove a certain threshold, in our case theta.Fortunately, logical expressions return 1 if theyare true, 0 if false. That is, (2>1) returns 1, (2<1) returns 0, or moreinterestingly (v>theta) returns 1 whenever v is bigger thantheta. This allows us to construct split functionsby multiplying the logical expression by its corresponding value and addingall such values together.

Arrays

SAY SOMETHING HERE???

Controlling XPPAUT from the .ODE File

Lines beginning with @ configure XPPAUT'sinitial settings. xlo, xhi, ylo, and yhi are probably the most useful settings as they let youset the initial window. dt requests a step size;the default is .05. totalsets the length of time to integrate over; the default is 20. SAY MORE???

Exercise: Implement a Hodgkin-Huxleyneuron with a constant applied current term as a parameter. Initially,m = m(Vrest),h = h(Vrest), andn = n(Vrest).With an applied current of 0, show that your neuron will not fire ifit starts with a membrane potential near rest but that it will fire if theinitial membrane potential is sufficiently above or below rest. Next, startwith the membrane potential at rest and increase the applied current. Whathappens?

Delay Equations

Synapses and other biological processes often act after a delay. Adetailed model might explain this delay by some set of chemical reactions,but often only the delay is relevant. XPPAUTlets us incorporate a delay with the function delay. For example, delay(x,3) returns the value of x at a time of 3 before thecurrent time. You also must set the @delay settingto an amount bigger than your longest delay.

In sample3, for example, we represented the delaydue to the build up and decay of neurotransmitters in the synapses by using adifferential equation for the synapse. Unfortunately, this gives us adifferential equation for each synapse. As a cruder approximation, we caninstead model a synapse as either on or off, switching after a delay of tau.This is sample4.ode:

# sample4.ode
# Post Inhibitory rebound in mutually coupled
# Fitzhugh-Nagumo cells; delay equation version.

# Parameter values. It's always better to give constants
# as parameters so you can change them without exiting
# XPP.
p vSpike=-1, eps=.2, gsyn=.5, theta=-.2,scale=10
p wshift=3.8, vshift=0, scale2=10
p tau=2, vsyn=-2

# Usually we talk of the FitzHugh-Nagumo model with
# f and g somewhat arbitrary, but to do numerics we
# need to pick specific functions.
f(a,b)=-b-15*(a-1)*a*(a+1)
g(a,b)=-b+scale*tanh(a*scale2)+wshift

# the slow system (equation) and the synapse obey
# the same rules for both neurons, so we use this array
# notation. The first line denotes the range of j (always
# j) and v[j] denotes the j'th v value, etc...
%[1..2]
w[j]'=eps*g(v[j],w[j])
s[j]=delay(v[j],tau)>theta
aux syn[j]=s[j]
%

# cell 1 receives input from 2 and 2 from 1, so we write
# these seperately since they are different
v1'=f(v1,w1)-gsyn*(v1-vsyn)*s2
v2'=f(v2,w2)-gsyn*(v2-vsyn)*s1

# initial conditions. Since I don't specify s values,
# they start at 0.
i v1=0, w1=-5, v2=-1, w2=-2

# configure XPP's options.
# delay=10 means that we must use a delay of less than 10
@ total=100, xlo=0, xhi=100, ylo=-1.5, yhi=1.5
@ delay=10
@ nplot=2, xp[1..2]=t, yp[j]=v[j]
done

Graphing it in XPPAUT using the defaultsettings gives us:

Note that the inhibition of a neuron happens a time of tau=2 after the other one begins to fire and is releasedtau=2 after the membrane potential of the otherneuron drops.

Auxillary Functions

Functions, in the sense described previously,cannot be graphed directly. Instead, we must create an auxillary function, asin the following:

aux syn[j]=s[j]

Recall that in our example, s[j] was a function.This line lets us plot synaptic levels just as we could (but did not) when welooked at sample3. Here, I plot both v1 (white) and syn1 (red) asfunctions of time.

Phase Planes

Let us return to our second ODE:

x1' = x2
x2' = −x1
x1(0) = 0
x2(0) = 1

We now plot X2 vs X1.

Use the Mouse to Select InitialConditions

Select Initialconds then either (M)ouse to select one initial condition by clickingwith the mouse or m(I)ce to repeatedly selectinitial conditions and plot trajectories. SAY MORE???

Flow

Rather than click over and over again, you can select Dir.field/flow and then (F)low to see the trajectories starting at evenlyspaced grid points in phase space. (You are prompted for grid spacing at thetop of the main window; simply press enter to select the default values.) SAYMORE???

Nullclines

Plot Nullclines by selecting Nullcline andthen (N)ew. Remember fixed points of atwo-dimensional system are precisely those points where the nullclinesintersect. SAY MORE???

Fixed Points

Points where the derivatives are all 0 arecalled fixed points; they may be stable, unstable, or nodes. To locate thefixed points exactly, use Sing pts. In oursimple example, we can just click (G)o afterthat. Select Yes when asked if you want toprint the eigenvalues (they will appear in your terminal window). A newwindow then appears with information about the fixed point it found. In ourcase, there is only one, but in general you may have to use one of the otheroptions in the Equilibria menu to find otherfixed points.

SAY MORE???

Exercise. Open sample3 in XPPAUT. Setgsyn to 0 to decouple theneurons. Select Graphic stuff then (D)elete last to go back to just drawing one plot(and hence one neuron's output). Graph W1 vs V1. What happens as you vary wshift? How does that affect the nullclines? The flow?Can you describe the role of eps?

Bifurcation Diagrams

The "AUT" in XPPAUT refers to AUTO, a program for drawing bifurcation diagramsfor ODEs which is included with XPP. SAYMORE???

Last updated on 10 July 2007.

Robert A McDougal: XPP Tutorial (2024)
Top Articles
Latest Posts
Recommended Articles
Article information

Author: Twana Towne Ret

Last Updated:

Views: 5964

Rating: 4.3 / 5 (64 voted)

Reviews: 87% of readers found this page helpful

Author information

Name: Twana Towne Ret

Birthday: 1994-03-19

Address: Apt. 990 97439 Corwin Motorway, Port Eliseoburgh, NM 99144-2618

Phone: +5958753152963

Job: National Specialist

Hobby: Kayaking, Photography, Skydiving, Embroidery, Leather crafting, Orienteering, Cooking

Introduction: My name is Twana Towne Ret, I am a famous, talented, joyous, perfect, powerful, inquisitive, lovely person who loves writing and wants to share my knowledge and understanding with you.