DSSRT: Dominant-Scale System Reduction Tool

Known compatibility issue: Pre-R13 versions of MATLAB do not support the predefined boolean constants used in this program. Make the following definitions at the top of each .m file: true = 1; false = 0;

Known bug: When 'P' command changes param file, time bar is not correctly restored and subsequent movement in time is not displayed with the vertical marker (nor are the markers).

Known bug: Occasionally, after switching windows, Matlab momentarily redisplays the DSSRT window with the time/position information incorrectly placed, and continues to update it in the correct place. However, the errant text remains on the window until the system is reset.

Possible bug: Lowering 'low-resolution multiple' parameter for AttEst too much makes computed "attractor" get smaller.

Known bug: If variable data .dat files are renamed, then reloading _FN.dat files referring to the previous names prevents the user from being able to change .dat files using 'O' command. The command will fail but DSSRT recovers.

Known bug: When functional net re-calculated, most recently calculated transition sequence still thinks it's "fresh" and command 'X' suggests its data is still being displayed in the main window, which it is not.

Version histories of individual components are listed in the associated .m files

2. Introduction to dominant-scale analysis

Nonlinear dynamical systems possessing hierarchical topological structures and multiple time scales arise commonly in the natural sciences. State-dependent coupling can add to the complexity, introducing new time scales that may not be explicit in the equations.

DSSRT embodies a computational technique known that can be used to reduce the study of such networks when expressed as a system of ordinary differential equations (ODEs), near a known orbit, to a set of low-dimensional approximate models. This method focuses on dominant scales in the relationships between variables during transient (far from equilibrium) dynamics. The technique adds rigor to intuitive reduction techniques that are ubiquitous in modeling high-dimensional coupled systems. In particular, it quantifies the robustness and parametric dependence of coherent temporal activity in more detail than a return-map stability analysis of periodic orbits, because it provides information along the entire length of the trajectory. The DSSRT software is somewhat tailored towards neural and chemical kinetics applications, although the underlying theory applies to a wider range of ODE systems.

DSSRT analyzes data that has already been numerically computed using a differential equation solver such as those in PyDSTool (http://pydstool.sourceforge.net), Bard Ermentrout's XPP software (http://www.pitt.edu/~phase), Neuron (http://neuron.duke.edu), GENESIS (http://www.genesis-sim.org/GENESIS), or any fixed time-step ODE solver from which you can export variable data to a file.

See the PDF documentation and the explanatory sections below, for more information about the theory behind the dominant-scale technique.

3. Getting started

3.1 Setting up the software on your computer

Make sure you comment and uncomment the appropriate lines in the file DSSRT_useroptions.m depending on whether you are running DSSRT under Unix/Linux or MS Windows, and change the directory settings according to where you've put your main DSSRT folder. Also, you can make sure that the `START_VERBOSE` option is set to `true` if you're new to DSSRT and wish to receive more system information in the DSSRT folder you are using.

In Matlab, change directory to the DSSRT folder. Type `DSSRT` at the prompt and enter the name of a directory containing network setup information (see Section 4. About the included demos regarding the included demos). Then type 'H' or '*' for further help within the main window that appears.

Variable data files (.dat) are assumed to have been generated at a fixed time-step in this version of DSSRT. You can generate these files using any fixed time-step ODE solver of your choice, provided you can output the data to a space separated text file in a column order that corresponds to that specified in a network setup (.cfg) file for your system. See the later section for details on the specification format, and the included demos for more information. Files included with the DSSRT system provide examples for use with the XPP software (http://www.pitt.edu/~phase).

3.1.1 Compilation of MEX files (Matlab R13 only)

Compilation of certain .m files (those known to exhibit speed improvements on compilation) can be done automatically the first time DSSRT is run, if a compiler is available to your Matlab installation. If you know that your Matlab compile is configured then set FORCE_NOTCOMPILE = true in the DSSRT_useroptions.m file, otherwise DSSRT will continue without compiling certain functions. Compiling other .m files by hand will most likely not improve the speed of the program (and will make some parts of the system crash on execution). I have tested the execution speeds of all components to optimize them either for the Matlab performance enhancement mode when interpreted (using the JIT system), or for compilation to MEX files, whenever they perform substantially better in those forms. This was done on a Windows NT system. If there is a problem compiling the the code on your system, the compile script may ask you to switch off the compile option in order for you to continue. The loss of performance would not be great in that situation, as the majority of the DSSRT functions work faster uncompiled when using the interpreter's internal speed enhancement (the Just In Time compiler).

3.2 Running the program

Start the program by changing directory within Matlab to the installation directory of DSSRT. To study a system that has already been set up (in its own sub-folder of this directory), say 'E_single_cell', type DSSRT E_single_cell at the Matlab prompt. This should initialize the program and load the model-specific configuration. Alternatively, typing DSSRT at the prompt by itself will bring up a dialog box in which you can specify the name of a system's folder.

3.3 Display a functional network diagram

A functional network diagram (defined in Section 6.1 Internal, external variables, and observables) is the most fundamental thing that you need to have in memory in order to use DSSRT to analyze your system. In most of the included demos (see Section 4. About the included demos below for details), there is a pre-computed functional network available for a pre-computed orbit of that system. That orbit is usually a stable limit cycle. To find out if the required file is present, press 'E' in the main window. By default, the filename in the dialog box will appear as 'demo'. Click OK, and if it is present (as the file 'demo_FN.mat' in the system's folder) then it will be loaded. Such a file is present for the 'E_single_cell' demo, for instance.

Alternatively, if the basic parameters were set up correctly in the .par file for the model system (see Section 7. Understanding the .cfg and .par file specification formats), and if the start and end times are set correctly (using commands '1' and '2'), then a new functional network can be computed. This is done using command 'U', which uses whatever resampling time-step is specified (by command '4'), and the current dominant-scale threshold (set using command '3'). The dominant-scale threshold is explained in Section 6.4 Dominance scale threshold, and re-sampling is explained in Section 6.6 Down-sampling the data (nOut).

3.4 Interpreting the graphical output

Firstly, you should familiarize yourself with some of the background theory in Section 5. Overview of DSSRT Matlab functions, and possibly in the theory document (see Section 3.7 Other documentation) for the more mathematically minded.

3.4.1 The ODE model as a network

Your ODE model system is represented graphically in the main window as a network of interconnected nodes. The .cfg file in the system's directory specifies the layout of the network in the window, among other things. Before a functional network has been loaded or computed, only the nodes of the network will be shown. These contain the names of the associated variables in the system. Small vertical "slider bars" will also be shown nearby to each node. These will show the current value of that variable using a small horizontal bar. If there is more than one slider beside a node, the additional ones will show the values of hidden ("internal") variables associated with that node, that are not visible to the rest of the network.

3.4.2 Color-coded arrows

After computing or loading a functional network various red and green arrows should appear on the diagram, which connect colored circles around the variable names shown. The arrows and the colors give information about which variables are "dominant" in the dynamics in a neighborhood of the computed orbit, at the time shown above the diagram. See Section 6.5 Active vs. Potentiated (a.k.a. potentially active) variables for details of the meaning of these graphics. The time-line is shown above the digital time / position readout, with a blue cursor marking the current time. Raw variable data associated with the current functional network diagram is plotted in the top pane of the main window for any variable of the system, and can be changed using command '6'.

But in the simplest terms, at a given time t, a red arrow from a variable 'x' to a variable 'y' in the functional network diagram indicates that x is having a "significant effect" on the time-course of y at that moment. A green arrow indicates that x does not have a significant effect on y, but that x could do so if it had a different value. Because DSSRT is told what range x can take (often it is a [0,1]-valued "gating" variable) this is usually telling you that if x was more highly activated then it could become a dominant variable on the dynamics of y. In the absence of any arrow between variables, there is either no physical connection present or at this time x could not take any value in its known range that would cause it to have a significant effect on y. The meaning of "significant effect" is described in Sections 5. Overview of DSSRT Matlab functions, 6.4 Dominance scale threshold, and 6.5 Active vs. Potentiated (a.k.a. potentially active) variables, and in the theory documentation.

3.5 Basic navigation

To navigate through the diagram as time changes along the loaded orbit, use the 'N' and 'B' keys for one time-step movements forwards and backwards (at the re-sampled time-step used to compute the diagram data), or ',' and '.' for ten time-step movements at once. Pressing shift with the latter two commands (i.e. pressing '<' or '>') will jump the cursor to the start and end times of the current diagram data (the very ends of the time-line).

3.5.1 Time markers

Time markers are useful for several functions within DSSRT. Think of them as defining a "selection" of time on which you wish to perform an action, as you would do by selecting text in a document editor that you wish to spell-check. Two markers can be set at any one time: one for a left end-point, and one for a right end-point. To mark a left time end-point, navigate to that time and press 'J'. This positions a small red tick on the time line. Navigate to the end time of your chosen interval and press 'K'. You can set the right marker before you set the left marker, provided you don't try to set a left marker at a later time position than a right marker that you have already set (and vice versa). You can change the marker positions by re-pressing the appropriate key at a different time position in the diagram.

Once set, you can navigate directly to a left marker by pressing ';' and to a right marker by pressing ':' (i.e. SHIFT plus the semicolon key). Notice that when you are exactly at an end-point either [L] or [R] will appear at the end of the digital time/position readout.

Positions of the time-markers are saved when the current functional network diagram is saved using command 'W', and they are over-written when one is loaded using command 'E' (or when a functional network is recomputed).

3.5.2 Switching to analysis of term sizes

After a calculation or file load of functional network data, the Greek symbol Psi should appear in the lower left corner of the main diagram window. This indicates that the source of the diagrammatic information about variable-over-variable "dominance" was derived from analysis of the "influence strengths". This is the default situation, but an alternative analysis of dominance is available. Dominance can also be measured simply by comparing the sizes of the input terms to the RHS of a differential equation. This is analysis is accessed by the command '-', which should change the Psi to the word 'Term'. The diagram states for the term size-based analysis were calculated at the same time as those for the "influence strength"-based analysis, and so the diagram state should switch instantly.

3.5.3 Other key commands

A complete list of the key commands is accessible by pressing 'H' in the main window, or separately in a file in the DSSRT documentation directory.

3.6 Compute a transition sequence

The next step of analysis using DSSRT is to compute an event transition sequence ('TS') for this functional network. This is a sequence of diagram states, each of which depicts the set of dominant interactions between the variables (see the definitions in Section 5.2 Transition Sequences). Usually you will want to restrict this to events involving a particular variable of focus (typically a voltage variable for neural models). First, check the variable of focus using command '5'. Mark out a time interval within which to compute the TS, as described above or you can use one loaded in a demo using command 'E'. Often you will want the selected interval to be one whole cycle period of a limit cycle, but it can be any sub-interval of the orbit you wish.

Now compute the transition sequence by pressing 'C', and press OK to restrict the sequence to the focused variable. (In very large systems, un-restricted transition sequences are slower to compute and do not highlight features specific to the local dynamics.) When this is done you have an opportunity to name your TS and add two additional data to be associated with it. For now just use the default values (by pressing RETURN repeatedly). Small red dots will appear above the time-line, indicating the positions of events found within the marked time interval. You can navigate through these event times cyclically by pressing the SPACE bar repeatedly.

You are now ready to try out other commands that aid your analysis of the structure inherent in your ODE system. Most of them require that a TS has been calculated before they can run. Use the help command 'H' to list them. The major ones have separate sections in this document explaining their use, e.g. the regime determination and attractor estimation / domain of validity utilities.

Please let me know if you have any problems or suggestions.

3.7 Other documentation

The `full` PDF documentation referred to in this document refers to a theory paper, DomScaleTheory.pdf, and a slide presentation that includes information about the attractor estimation algorithm AttEst.m. Beware that some of the notation in those documents may not yet be entirely consistent with that adopted here.

4. About the included demos

4.1 Single excitatory Hodgkin-Huxley cell with no synaptic input

The `E_single_cell` demo is the most basic demonstration of the functional network (event, epoch, regime) analysis applied to a neural model. The cell regularly fires with a period just under 50ms due to a drive current. The regimes determined from this model agree with those derived by formal methods in a paper by Suckley and Biktashev (see the PDF theory document for details).

4.2 Excitatory and Inhibitory two-cell network

The demo 'EI_demo' consists of two hippocampal neurons E and I (modelling a pyramidal cell and an interneuron) coupled through chemical synapses. The EI.ode file for XPP contains the set-up for this network. The principal behaviour of this system for the parameter regime given is a synchronous state in the gamma frequency regime (30-100Hz), in which the I-cell fires twice per cycle (known as a `doublet`) compared to the E-cell, which undergoes a temporary suppression every half cycle due to the incoming strong inhibition. Various data sets of the variables produced by XPP are listed, as the I->E coupling variable g_ie is varied around its original value. The saved functional network demo_FN.mat contains one period of data for a functional network for one of those data sets.

4.3 Single excitatory cell with inhibitory synaptic input

The demo 'E_demo' is a stripped-down version of 'EI_demo', for a simpler view of a single neuron. In particular it more explicitly shows the roles of the action potential currents. The inhibitory cell has been replaced with a periodic inhibitory synaptic input. A sample variable data file, and a functional network with markers set for one complete period, is included.

4.4 Fitzhugh-Nagumo oscillator

'FHN_demo' takes a look at applying the dominant-scale technique to a non-neural equation. The FitzHugh-Nagumo equations are set up here in a classic relaxation oscillator regime. In this case, there is a 'synaptic'-like input added. In one form of the two FHN systems provided, the dummy `w` variable is dropped from the analysis (see the XPP .ode files for details). The AttEst function does not provide reliable information about FHN orbits at this time, because this example system involves non-dissipative dynamics, and AttEst does not yet deal with this class of ODE system.

4.5 Enzyme kinetics model

This demo is a simple Michaelis-Menten model, showing how an ODE system with more complex inputs can be re-arranged in order to be amenable to the configuration specifications of DSSRT. The dynamics are very simple however, and the demo is to illustrate that the correct regimes are resolved from the orbits, in correspondence with classical two time-scale analysis of this system. (Strict `speed fussiness` in the RegimeDet algorithm is required here.)

5. Overview of DSSRT Matlab functions

Please refer to the full documentation PDF file for more detail about the theoretical concepts discussed briefly below.

5.1 Functional Networks (FuncNet.m)

The Matlab function FuncNet.m is run using key command 'U'. It calculates subsets of the set of observable variables in the system for each time-step over a selected time interval. At each time step the subset contains the variables are most dominant in affecting the asymptotic target of the system. This is judged according to a `dominant-scale threshold` i.e. a scale of influence of one variable on another. The set of all the unique time-ordered subsets can be viewed as a sequence, defining the states of the functional network's diagram. The dominance of variables on others is sometimes non-intuitive. Our definition of dominance is given in Section 6.3 Influence strength (Psi).

Instead, we define a measure of the effective degree of influence that a change in an `input` variable has on the asymptotic `target` of the observable is ranked against all other forces exerted by the `input` variables to the RHS of the observable's equation, and the ranked set (in descending order) of variables having the largest such influence is truncated after their size becomes a sufficiently small fraction (the dominant-scale threshold) of the most influential variable.

More precisely: A type of sensitivity calculation determines the temporal pattern of the “dominant” variables over an orbit (e.g. an attractor), according to the definition of scale of influence. This partitions the phase space in the temporal direction into epochs, according to a user-defined dominant-scale threshold. The subsets of dominant variables in each epoch define the transition sequence. Within each epoch, the local basin of attraction to the invariant manifold is approximated by a lower dimensional manifold corresponding to only the locally dominant variables.

5.2 Transition Sequences

A Transition Sequence (abbrev. TS) is a set of states of a functional net that describe all changes in the network's dominant variables over some time interval, i.e. excluding the contiguous non-changing states that typically occur in the original functional network's state sequence. Each transition is called an event, and marks the beginning of a new epoch, during which a certain set of variables are `active` (most dominant), and those remaining are one of two types of inactive variable: potentially active, or non-potentially active. The first set is known as the Actives set, the second is the Potentials set, and the remaining variables are only implicitly defined by exclusion from the first two sets. (Additionally, the user has the option to distinguish changes in the order of the Active sets as indicating a transition, using command '7'.) The most recently computed TS is known as the 'fresh' TS. Any saved TS can be loaded in the current version of DSSRT, regardless of whether it was computed for the same parameters as currently set, or even for the same variable data set. Until this irritating matter is addressed in a future version of the program, a temporary feature is that many functions requiring a TS specify that the TS be fresh. This guarantees it to be consistent with the whole ODE system currently set up.

A set of related epochs make up a reduced dynamical regime, within which a certain set of equations for dynamical variables is most appropriate for analysis, according to standard multiple-scale analysis techniques (see Section 9. Domain of Validity / Attractor Estimation for automated features).

Different transition sequences can be calculated for different saved data sets that contain different time-evolutions of the ODE system's variables, e.g. arising from different initial conditions or parameter changes in the original ODE system. A notion of `distance` from a transition sequence (e.g. a periodic cycle) can be defined in terms of the sum of the errors in comparing a similar sequence of the same length to the original, according to a particular algorithm for making similar sequences the same length. This leads to speculation of the meaning of considering TS's in the same fashion as discretized dynamics, and thus of the stability and bifurcation of transition sequence cycles.

5.2.1 Viewing and comparing transition sequences

Command 'V' compares transition sequences of equal length, in terms of their active variables only (their potentially active states are ignored). If two TS's are not originally of equal length the command 'Z' will attempt to pad the shorter one in an intelligent way so that the positional errors between them are minimized. Command '#' controls the fussiness of the padding algorithm, although it is still in development.

Key to the horizontal `bar` display of viewed transition sequences (using command 'V'):

· Contiguous blue lines indicate original portions of the sequence.

· Broken lines in black indicate portions that were padded by the sequence matching algorithm in order to minimise the number of errors between sequence members over the length of the sequences, and to make the sequences the same length. (Note that if the transition sequences were focused on a particular subset of the variables, then only differences involving those variables will be indicated.)

· Remaining errors in sequences generated by the matching algorithm (see StepNet) for a sequence are shown by red markers at their sequence positions.

5.3 Main window interface and commands (StepNet.m)

StepNet.m provides the principal user interface to control the derivation and comparison of functional network diagrams and their transition sequences.

5.3.1 Information at the current time position

Command ‘`’ lists information about the variables at the current time. It includes their values, their time scales of reaction, and some rough indicators of the accuracy of the reduced set of variables indicated in the diagram (these are only experimental at this time). The first indicator uses the following identity that applies to conditionally linear systems:

where R is a remainder (error) term. This equation is simply the current balance equation for V applied to the value of V at V_\infinity. The value |R|/(N-|Acts|) measures the average magnitude of each input term’s contribution to the remainder.

Command ‘@’ lists information about the dominant scales at the current time position, ranking the variables according to their influence strengths. Information is also provided about the potential influence strength rankings and values.

5.4 Specification of differential equations

Differential equation specification is achieved using the DEQNS, GAMMA1, GAMMA2, CFAC and DEPARS specs in the .cfg setup file.

Take the example of a Hodgkin-Huxley neuron, such as that in the included demo `E_demo`. Typically, a quantity such as `mtau_recip` (governing the relaxation time of the m variable) is not a constant value multiplying the gating variable m in a right-hand side, but a self-contained activation function. Thus, when they are identified as such functions, argument 3 = `u1` indicates what `mtau_recip` is a function of (note only 1 argument is allowed for those functions). In such cases the power in argument 4 is ignored (and so typically set to 1). `mtau_recip`, etc. are known to DSSRT as functions of other variables, and not just constant parameters, because the INPUTS specification for the associated gating variables are non-empty. (DSSRT takes these names and uses the MATLAB interpreter eval to execute functions with these names, and so .m files having those same names must reside in the ODE system's folder.)

CFAC specifications are optional. For a differential equation, this specifies the name of a DEPAR (a constant) that multiplies the derivative. For neural equations this would be a capacitance. If the spec is omitted, the value defaults to 1 (this is in units of micro Farads / cm ^ 2 for neural equations).

For each differential equation, there should be no more than INPUTS(eqn) number of GAM1TERMs + GAM2TERMs. Also, for fewer than that number, those present in INPUTS that are not in a Gamma set may only be internal variables. See the examples provided in the demos.

5.5 Reduced regime determination (RegimeDet.m)

This function takes an event transition sequence and consolidates the epochs into a reduced set of regimes, which are more amenable to analysis. The function does this by looking for commonality between epochs, according to a simple set of rules. It also explicitly identifies fast and slow variables with respect to the focused variable, and predicts appropriate quasi-static bifurcation parameters for each regime. Section 8. Reduced Dynamical Regimes contains more information.

5.6 Domain of validity / attractor estimation (AttEst.m)

Although still in development, AttEst.m uses an efficient algorithm to calculate the domain of validity for a transition sequence or sequence of regimes. Using slightly different settings, it computes a conservative estimate of a cross-section of the basin of attraction around a known orbit, so it can be used to study the orbit's robustness to perturbations. This can be done for any of the external variables. More information is given in Section 9. Domain of Validity / Attractor Estimation.

6. Definitions and descriptions of main parameters

Please refer to the full documentation PDF file for more detail about the theoretical concepts introduced below.

6.1 Internal, external variables, and observables

These are subjective definitions which depend on the interpretation of the structural components of the system. On the RHS of an ODE for a certain variable `V`, if certain terms are considered due to an input `external` to the sub-system that V represents, then the variable contributing to that input is an external variable. Conversely, if inputs to V arise from dynamics that are considered internal to V's subsystem, then that input's variable is an internal variable. E.g. for the Hodgkin-Huxley equations, we might consider the sodium and potassium spiking currents to be internal to the cell, and so variables m, n, and h are internal variables, whereas a synaptic input from another cell means that it's variable s is an external variable. If we consider the equation for the synaptic variable existing in a separate sub-system to the soma then the pre-synaptic membrane potential is an external variable since it inputs to the equation for s. The synaptic equation has no internal variables in the absence of adaptive processes (e.g. depression, facilitation). An `observable` is merely a synonym for an external -- it is a variable that can be `observed` by other sub-systems.

6.2 Passive and dynamic variables

We call variables that have an associated differential equation dynamic variables, and those that do not are passive (non-dynamic). Passive variables play the same role as non-passive variables in input terms to a differential equation, but they might be truly formal (i.e. actually non-varying) or be an externally-varied signal. This terminology is used to distinguish this issue from that of dominance, where the term `active` is used (see Section 6.5 Active vs. Potentiated (a.k.a. potentially active) variables).

6.3 Influence strength (Psi)

Loosely speaking, the influence strength of an input to some variable V's ODE is the amount of control that the input has over the asymptotic target of the variable. The usefulness of this quantity in DSSRT relies upon certain assumptions about the ODE for V, as discussed in the full documentation (PDF file). In particular, that trajectories V(t) contract sufficiently quickly towards the asymptotic target, that the motion of the target is a reasonable approximation to the actual trajectory, in some specific sense.

6.4 Dominance scale threshold

This is the critical value defining how exclusive or inclusive you wish your sets of active inputs to each observable variable to be. It is set using command '3' (which also sets the time scale threshold discussed in Section 6.7 Time scale threshold, `fast` and `slow` variables). Mathematically it defines the `Order 1` scale of variables in the system, and is the reciprocal of the ubiquitous small parameter epsilon used in multiple-scale analysis. Formally, epsilon is the fraction of the largest ranked influence strength (for any observable), such that any other input's influence strength less than this critical fraction will not be included in the set of active inputs acting on a the observable variable, at one moment in time (see the full documentation for a mathematical definition).

Increasing epsilon away from 1.0 will loosen the tolerance on the strength of influence of other candidate-active inputs to the observable, and allow weaker-influencing inputs to become considered `active'. Decreasing towards 1.0 will tighten the tolerance on influence strength and ignore more and more of the weaker-influencing inputs. Values of 1 -> 1/10 give different intuitive perspectives on what is `dominant` at any given time, at least for ODEs describing neural membrane channel dynamics. Plotting the influence strengths along the orbit loaded in DSSRT (using command '~') indicates where to look for "spectral gaps" that are appropriate places to set your epsilon value. Systems with many variables may have multiple gaps, and these tend to highlight the fact the system may be driven by some variables, but modulated by others.

Note that in the accompanying PDF documentation, the epsilon scale threshold is currently referred to only by its reciprocal, and then it has the symbol epsilon.

6.5 Active vs. Potentiated (a.k.a. potentially active) variables

Green-coloured links in the diagram (and certain special colours for a subsystem's internal variables) show when an input has the potential to be in the set of most dominant variables acting on some observable, if it is not currently in the set of actives. This can only be defined using knowledge (or estimation) of the maximum absolute value of that input, and is somewhat sensitive to the actual choice of that value (unlike for the actives set itself, which is insensitive to this choice provided it is sufficiently conservative).

The procedure for determining this at any time step is to find the maximum value of the influence strength for the variable's gating coefficient s_i over the range 0 -> 1, given the actual current values of all other variables. This value is then compared to the current set of actives, and if it would belong to that set (using the discrimination of the scale threshold), then that input is considered potentiated. This information is less useful than knowledge of the actives set, since it is only a guide to the potential influence that variable could be having at some time, i.e. if it were to exactly reach it's point of maximum influence. Gating coefficients with associated cross-multiplying variables are not varied as part of the maximization, even though in reality changes in the primary variable might always be coupled with changes in the cross-multiplying variables. This is another reason why the set of potentially active variables should be used merely as a heuristic guide to where to make more detailed analyses.

One helpful trick with making the most from 'potentials' is to do the following. If you are studying behaviour near a known orbit, observe the range of variable values taken on that orbit and in perturbations from it that you might be interested in. Then, restrict the BOUNDS for that variable accordingly in the .cfg specification (although leaving a small safety margin is recommended). That way, the accuracy of the green 'potentiated' arrows is improved, having been made more relevant to this dynamical regime. Be warned however that significant changes in system parameters may invalidate these bounds and you will have to make sure you adjust them accordingly.

6.6 Down-sampling the data (nOut)

If the sample resolution of the variable data in the .dat file was set too small during the original numerical integration of the ODEs, DSSRT may take a long time to load files and perform major calculations. Use the nOut setting (command '4') to re-sample data every nOut steps when calculating the functional network to improve speed and efficiency. However, beware that overzealous down sampling can lead to quantization artifacts, namely missed or lumped transitions in the state sequence. Such artifacts can make comparison with other data sets unreliable. The visible resolution of the variable displayed in the top pane of the main window is also determined by this parameter, and can be used as a rough indication of whether it is appropriately set for a particular dataset.

6.7 Time scale threshold, `fast` and `slow` variables

This is also set using command '3', along with the dominance scale threshold. It is only used by the attractor estimation algorithm (see Section 9. Domain of Validity / Attractor Estimation) and for generating reduced dynamical regimes (see Section 8. Reduced Dynamical Regimes). In the PDF documentation this threshold is referred to only in its reciprocal, using the symbol gamma.

The threshold is used to classify a variable within a particular epoch as either possessing an "Order one" speed of response to perturbations, compared to the focused variable of an analysis (e.g. the membrane potential of a neuron), or a `fast` or `slow` response compared to the focused variable (order gamma or order 1/gamma, respectively). These classifications are the basis for further reducing the number of effective dynamical variables in local models.

Note that we do not consider variables whose time-derivatives are very small or very large to be necessarily slow or fast. For instance, along or near an orbit such as a limit cycle, variables with very fast responses can be slaved to a slower variable, so that their time-derivatives remain small, but their response to perturbations would be very fast. DSSRT uses an additional classification -- `near-constant` -- to describe any dynamic variable that happens to not significantly vary over an epoch (relative to the variable bounds declared in the .cfg file).

7. Understanding the .cfg and .par file specification formats

7.1 Overview

· A .cfg file is a list of specifications for an ODE system and its display in the DSSRT window.

· Comment lines begin with `#`. All other non-blank lines will be treated as specification commands.

· Command lines must not end with comments or additional text.

· Variables not generated by explicit DEs must be declared internal variables (e.g. `u2` variable in `E_demo`), e.g. use for 'external' time-varying inputs in non-autonomous systems. Strictly, such formal 'internal' variables are still 'observables' (see definitions in Section 6.1 Internal, external variables, and observables).

· The order of declared variables must be the same as the order of the variable data file columns.

· .dat data files may contain more columns than the total number of system variables, for instance if the software that computed them also stores auxiliary data in those columns (e.g. in XPP). This is true for some of the data files in the included demos. Any additional columns are simply ignored by DSSRT.

· A DSSRT diagram state is defined to be the states of the set of nodes and the links between them.

· Non-regular state maps (see Section 7.4 Additional node states and state maps below) must be used to include activity of internal variables that are not given explicit nodes (e.g. the action-potential gating variables in `EI_demo`).

· The DSSRT diagram drawing area has normalized co-ordinates [0,1] x [0,1], for use in the diagram specifications.

· A .par file is an unlabelled list of values for various basic parameter defaults loaded at start-up. The order of the entries in the file are as follows:

· CFG file name, var data filename, start time, stop time, dominant scale threshold, variable data resampling rate, ignore actives/potentials order switch, period guess (if applicable), margin % error of period in searching, time scale threshold

The configuration file for a system can be created semi-automatically by running the DSSRT_utils.py utility in PyDSTool’s Toolbox directory, on a version of the system expressed in the syntax of PyDSTool. (I am happy to help create your configuration file for you.) The rest of this section deals with the structure of the .cfg file format.

7.2 Correct order of system specifications

(1) Variable declarations

VARSEXT

VARSINT

The order of the variables listed { VARSEXT, VARSINT } must correspond to the order of numerically integrated variable data in a .dat file loaded into DSSRT.

(2) Network connectivity, i.e. variable inter-dependence

INPUTS

One for each variable, in any order. A differential equation must exist for any observable having more than one external variable input.

(3) Information about variable bounds

BOUNDS

UNITBOUNDS

There can be a number of these specifications, and neither is required, provided all variables either have their own BOUNDS spec or are mentioned in the UNITBOUNDS list. Due to the common occurrence of unit-interval variables, this list obviates the need for individual BOUNDS specs for each of these variables.

(4) Differential equation set-up

(4.1) DEQNS

One of these specifications, listing which variables have an associated differential equation being declared to DSSRT.

(4.2) GAM1TERM

GAM2TERM

At least one of either, for every declared differential equation. These declare all the terms in the right-hand side of the differential equations, according to the classification of Gamma1 and Gamma2 term-types (see the theory paper for details). `Tau-reciprocal` and asymptotic `target` values for each term can either be declared DE parameters (see DEPAR specification), or .m files in the set-up directory that specify the appropriate functions to generate these values (see Section 7.6 User-defined .m functions for use in DE specification). Additionally, Gamma1 `targets` can be the names of other external variables of the system. (For instance, this permits the modelling of electrical synapses in neural equations.) Numerical values for the `tau-reciprocal` and `target` values cannot be used directly in a GAM1TERM or GAM2TERM declaration, and no arithmetic combinations of parameters are possible in the declaration (ie. separate DEPARs for every coefficient are needed even if they are algebraic combinations of the fundamental system parameters!). Inputs in the Gamma1 set are allowed no more than one cross-multiplying internal variable (e.g. the 'h' inactivation variable in the sodium activation term m^3 * h, for the Hodgkin-Huxley equations). This should be enough for most modelling purposes. Currently, Gamma2 terms cannot include cross-multiplying variables, but this will be supported in a later version.

(4.3) CFAC

Optional constant multiplicative factor on the derivative of the DE named in the first argument. Defaults to 1.0.

(4.4) DEPARS

As many of these as there are parameters used in the GAM1TERM and GAM2TERM specs

(5) Functional network diagram set-up

(5.1) NODE

NODELABEL

NODEADDSTATES

NODESTATEMAP_ACTS

NODESTATEMAP_POTS

One set of these specifications, in this order, for each external variable. If no additional states for the node are needed (the usual case) then these specifications must still be present, using `NODEADDSTATES <Var> 0` and then `NODESTATEMAP_ACTS <Var>', `NODESTATEMAP_POTS <Var>', where <Var> is the relevant observable variable name. See `EI_demo` for a typical use of the additional node states, in order to avoid the clutter of too many nodes for the action-potential generating variables in the Hodgkin-Huxley equations. Note that a node label does not have to correspond to the variable name used internally to DSSRT.

(5.2) LINK

One of these for each variable interdependency that involves variables with their own diagram NODEs (i.e. for every coupling between observable variables described by terms in a differential equation). This specification displays an arrow between the connected variables in the main graphical display, which will be colour-coded (or absent) depending on the 'functional network' state (see Section 7.5 Links). The display of unwanted arrows can be suppressed by specifying `1` in the optional 7th argument. These specifications can be in any order. (A future version will alleviate the need to specify explicit arrow co-ordinates, although these can be generated automatically using a Python code that I have recently developed.)

(5.3) VBAR

One of these for every declared variable, in any order. The final argument is the switch for allowing log-scaling: 0 = `off`; 1,2 = `on` ... but 1 means the low-value end is magnified, and 2 means the high-value end is magnified (e.g. for variables that usually rest near their maximum values rather than their minimum, such as the sodium inactivation variable `h` in the Hodgkin-Huxley equations).

7.3 Full syntax for DSSRT specifications

VARSINT <varname> [ <varname> ... ]

VARSEXT <varname> [ <varname> ... ]

INPUTS <dest_var> [ <varname> ... ]

BOUNDS <varname> <lo_val> <hi_val>

UNITBOUNDS <varname> [ <varname> ... ]

DEQNS <eqn_subject_name> [ <eqn_subject_name> ... ]

GAM1TERM <eqn_subject_name> <taureciprocal> <extvar> <power> <target> [ <intvar> <power> ]

GAM2TERM <eqn_subject_name> <taureciprocal> <extvar> <power>

CFAC <eqn_subject_name> <parname>

DEPAR <parname> <value>

NODE <obslabel> <x-centre> <y-centre> <radius>

NODELABEL <obslabel> <name> <x-textpos> <y-textpos> <textsize>

NODEADDSTATES <obslabel> <numstates> <style> [ <style> ... ]

NODESTATEMAP_ACTS <obslabel> [ <statemapfunction> ... ]

NODESTATEMAP_POTS <obslabel> [ <statemapfunction> ... ]

LINK <source_obslabel> <dest_obslabel> <x-tail> <y-tail> <x-head> <y-head> [ <visibility-switch> ]

VBAR <varlabel> <x-pos> <y-botpos> <y-height> <variablelo_val> <variablehi_val> <log-scale-switch>

Notes: <obslabel> = a declared external (observable) variable name

<eqn_subject_name> = a declared variable name

7.4 Additional node states and state maps

Additional node states can be specified using NODEADDSTATES. These states have values corresponding to their position in the command arguments. DSSRT always includes an obligatory state 0, meaning all inputs are inactive and unpotentiated. In the `EI_demo` example (where the sodium (Na+) gating variable is mu1, the potassium (K+) gating variable is nu1):

NODEADDSTATES u1 6 g- r- y- k- r: y:

'g-' = either Na or K potentiated (state 1), 'r-' = Na active only (state 2)

'y-' = K active only (state 3), 'k-' = both currents active (state 4),

'r:' = Na active, K potentiated (state 5), 'y:' = K active, Na potentiated (state 6)

These codes are standard Matlab plot styles: a letter for the colour name, followed by a token giving the line style. The state meanings all refer to states of internal variables associated with the current observable (u1) only, because these variables are not separated out to have their own links to u1 showing those states individually in the diagram.

The possible transitions of these states are defined in the node state maps, of which there are two in this example, since there are two internal variables in the candidate actives list for u1 implied by INPUTS (the `m` and the `n` variables mu1 and nu1). The state maps must have 2 columns (previous state on left -> new state on right), and same number of rows as total number of states (although some will be only for potential states, or for active states). Unused entries in the domain of the maps should map to `-1` by convention. The maps are names of Matlab m-files in the same directory as this .cfg file. These also have to be specified in the same order that the associated internal variables are listed in the VARSINT specification. There is no other way for DSSRT to know which map is associated with each variable.

NODESTATEMAP_ACTS u1 nodeStateMapActs_M nodeStateMapActs_N

NODESTATEMAP_POTS u1 nodeStateMapPots_M nodeStateMapPots_N

7.5 Links

'Links' are the connections between the nodes of the network, i.e. the couplings between observable variables (these cannot include internal variables). The start and end points of arrows appearing in a functional network diagram are given by of LINK specifications. Each LINK specification should correspond to an external variable input (arg 1) to an observable (2nd argument), as defined by a relationship in an INPUTS command. All links have default state maps, supplied by DSSRT automatically. The four arguments following are the x1, y2, x2, y2 values of the arrow start & end, respectively.

There is an optional final argument in all these LINK commands: if it is 0 or not present then the link is displayed, but 1 makes the link invisible. This is useful when there are links that you don't care about and you don't want to clutter the diagram. DSSRT dictates that every interaction between two observables (specified by INPUTS) must have an associated LINK command, and so for the unwanted links, specify the 1st and 2nd arguments as usual, but specify any co-ordinates (for instance, zeros) for the arrow and put a 1 at the end!

7.6 User-defined .m functions for use in DE specification

DSSRT looks for user-defined .m functions in the system's set-up folder if the names specified in a GAM1TERM or GAM2TERM specification for a differential equations are not known as parameter names (in a DEPAR specification). Currently, these functions can only be functions of the single variable declared as the variable name <extvar> (see the full syntax specifications). This restricts the type of equations that can be declared in DSSRT at this time.

You can adapt the .m files for the asymptotic (x_\infinity) function xinf.m, and the time scale function for x (specified as a rate, i.e. reciprocated), in xtau_recip.m, where `x` is `m`, `h`, etc. in any of the neural demos. An annoying feature of some of the standard functions found in the neural modeling literature is that some constituent activation functions for xinf and xtau_recip have singularities in their usual range (e.g. m_alpha(V), called `ma` in the function minf(V), for E_single_cell). Care must be taken to add a special case in such functions to catch any possible division by zero. So watch out for those...

Also note that you don't have to define your xinf and xtau_recip functions in terms of forwards and backwards rates. Specify them however you wish. They only have to be self-contained functions (all other parameters and states in DSSRT are not observable by them).

8. Reduced Dynamical Regimes

8.1 Overview

Standard multiple-scale asymptotic analysis can determine `regimes` within which a reduced set of equations is appropriate for accurate study of the dynamics, according to various constraints (and in particular, over restricted time intervals). This analysis generally makes use of the time scales associated with each variable, and involves a lot of intuition and artistry to get right. The algorithm incorporated in the Matlab function in RegimeDet.m aims to encode some of the intuitive and analytical steps necessary to help you quickly determine a set of reduced dynamical regimes for your system of ODEs. The algorithm brings together information from several neighbouring epochs and establishes whether they should be considered together as a single coherent regime, according to a few very simple rules about the types of event at epoch transitions. RegimeDet identifies fast time scales along an orbit and then looks for slow variables that can be used as quasi-static parameters for the purposes of bifurcation analysis within the regimes. While the code is designed to make these choices intelligently and robustly, it has only been tested on a few example systems and may not always give the most appropriate results. At this time, the code also only applies to regimes associated with a single variable (see the published journal paper for details). Please contact me with comments and questions so that I can try to improve the functionality of this algorithm.

8.2 Using the RegimeDet function

Reduced dynamical regimes can be computed from a fresh transition sequence using command 'R'. Keeping the `ignore order of actives` switch off (command '7') may help the algorithm to detect time scale changes more accurately. Set the verbose option on (command '?') before calculating the regimes in order to see the detailed analysis of each epoch before the final regime information is displayed. Use this option to help guide your fine tuning of the time scale threshold (command '3') if there is "chatter" between epochs creating unwanted extra regimes. As always, your intuition will be an equally important tool in getting sensible results.

It is common to be studying a limit cycle, whereby the transition sequence used for the regime determination is itself cyclic. Specifying this when asked will help the program to identify the correct starting place of regimes in the cycle, regardless of where your TS begins and ends. (Thus the first regime's displayed time interval may not correspond with the beginning time of your TS.) For non-cyclic data, no information is available for epochs before the first one in your TS, or after the last one. In that case the system will begin a new regime with the first epoch, and will not be able to suggest a quasi-static bifurcation parameter for the final regime.

In the report about the regimes that is displayed in the console, parentheses next to variable names indicate whether those variables have been deemed fast (`F`) or slow (`S`) in their response, relative to the principal focused variable, according to the time scale threshold. Any cross-multiplying variables associated with observables in an input term of a differential equation are marked with `X`. These may or may not be internal variables depending on your set-up. Fast variables are slaved to their input, slow variables are effectively constant over the regime. What constant that is depends self-consistently on the motion of relevant variables in previous regimes, but is not suggested automatically in this version. The last value of a variable from a previous epoch on the original orbit could be used, before it went `slow`. Alternatively, a value from the middle of the present epoch, again lying on the original orbit, could be used.

Variables that are neither fast nor slow in their response may nonetheless change by only a tiny amount over an epoch or regime. By default, RegimeDet will track this according to a "small variation threshold" (initially set at 5% of a variable's declared bounds). If a variable does not change by more than this fraction of its bounds over an entire regime then it is classified as `near-constant`, and is marked with a `C` in parentheses in the report. This is in addition to any other classification that variable receives in that regime. Near-constant variables cannot be used as quasi-static bifurcation parameters for that regime if they are also classified `slow`, because they are presumed not to change enough in the regime to cause a bifurcation. This presumption must always be verified by the user, because it is not always valid. The threshold can always be lowered to zero if this classification is not desired. Because this classification is not measured relative to any other variables, DSSRT includes the focused variable in its checks for being `near-constant`. However, it will not classify this variable as such in the final regime report, because it is assumed that the focused variable is the usual variable in which to consider perturbations. Instead, DSSRT mentions the status in the verbose report made during its calculations.

In formal analysis of a regime the classifications `fast`, `slow` and `near-constant` may be violated when studying nearby orbits. This should be taken into consideration when setting up the most suitable model for each regime. DSSRT will provide some guidance about the domain of validity for a chosen model with respect to perturbations, etc., if the AttEst algorithm is used (see the next Section). Remember that the output of the RegimeDet function is only meant to be a guide for your analysis: not all of RegimeDet's suggestions can be used together. For instance, a variable may be noted as changing less than the `small change threshold` throughout a regime, but its timescale may not be slow.

8.3 More information about the algorithm

8.3.1 When does a new regime get started?

The algorithm's use of the internal `speedChange.flag` assumes that only one variable will change timescale status at a time, otherwise only the last variable to be checked will be correctly updated. Therefore, a small enough time-resolution of variable data is assumed.

An event in which one or more dynamic variables leave the active set may cause a new regime to be started with the new epoch after the event. This will not be done if there are no dynamic or quasi-static variables present in the regime being ended. The idea behind this is that a regime should include a means to move on to the next regime built into it (assuming that we are studying non-constant trajectories!). An overly-reduced model may contain a fixed-point that does not move slowly (quasi-statically) over a longer time-period, in order that some bifurcation or other transition occur in order to continue along the original orbit. In the absence of dynamic variables or quasi-static parameters, this failure of the reduced model is most likely.

8.3.2 Algorithm options: `fussiness`

Changes in timescale status in between events may go unnoticed without the `speed fussiness` option set. However, this might be fine for your problem, depending on the stiffness of your equations. For coupled neural systems, it seems so far that you needn't be `fussy` very often.

Note that if the speedFussy option is set at `not strict`, the algorithm will flag the start of a new regime if a variable changes its status between `near-constant` and not `near-constant` between epochs. If you set the option to be `strict`, the algorithm additionally does not allow time derivatives for any variable larger than a scaling of the specified threshold. The scaling is proportional to the specified bounds on the variable.

There is a `fussiness` option for what we could call `fast leavers` from the actives set. Normally the algorithm will start a new regime when a non-passive (`dynamic`) variable leaves the active set, provided there are either dynamic variables or quasi-static bifurcation parameters in the regime already (see section 8.3.1 When does a new regime get started?). With the `fast leavers` option set, a new regime will only be started if a variable that is leaving the actives set is not fast. The rationale for this is that keeping a fast variable in a regime does not increase its dynamic dimension, because the dynamics of that variable would be approximated by its asymptotic behaviour (an `adiabatic elimination`). Therefore this helps to simplify your regimes, and create fewer of them at little additional expense. However, it is sometimes the case that you really want those fast variables to go away when your intuition tells you that they are no longer needed, for instance: when they would be being kept over very long epochs which otherwise would have much simpler reduced dynamics. Both options should be tried, and the original transition sequence checked over to see what seems most appropriate. Neither option is necessarily more correct.

As an example of the `fast leavers fussiness`, more meaningful regimes are generated for the uncoupled Hodgkin-Huxley cell in the `E_single_cell` demo with this option off, whereas the opposite is the case for the `Stellate_single_cell` demo. This has partly to do with there being more variables in the latter demo, and so it is more appropriate to try to reduce the number of regime changes.

8.3.3 Slowness

Variables that are initially deemed slow in their response for a whole regime are further checked: their (instantaneous) characteristic timescale values at every time point are checked against the growing duration of the regime. If they are of the same order, i.e. according to the time scale threshold, then the variable is upgraded to a `normal speed` variable. This occurrence is reported in the verbose output of the algorithm. Use of this will depend on the stiffness of your system.

8.3.4 Dimension of reduced regimes

There are three versions of a reduced regime's dynamical dimension calculated. The first (A) is the number of intrinsically dynamic variables active in the regime i.e. those that have associated differential equations. The second (B) is the number of dynamic variables minus those variables that have been classified `fast` or `slow`. The third (C) is the same as B but also subtracts out the variables that remain `near-constant` throughout the regime.

8.3.5 Regime transitions, bifurcations

In future versions of DSSRT, basic fixed point stability analysis will be automated from the regime information, and it might be possible to generate a "pertinent" bifurcation sequence for the regimes. The algorithm will currently not tell you whether a regime is over-reduced, in that it might have a fixed-point with no quasi-static variables to cause that fixed-point to ever bifurcate or for a new regime to be otherwise entered.

8.3.6 Practical hints

Finally, please remember that this algorithm is very much in development. Try adjusting your parameters and options when using it, and see if you can find a set of regimes that make sense to you. If you get stuck or find annoying issues crop up, look at the code and/or contact me. I will be very happy to investigate further, so that the algorithm can be made more robust.

8.4 Phase-plane diagrams and movies

For many physical systems possessing multiple scales and modular structure, it is common for a small number of variables to be the effective dynamic variables at any point in time. At the end of the regime determination, you are given the option to display a sequence of phase-plane diagrams that characterise the dynamics through the regimes. The sequence of diagrams might elucidate bifurcations occurring in a regime, for instance.

If the effective dimension of a regime does not happen to be 2, you are presented with a choice for the y-axis variable in that regime. (The x-axis is automatically selected to be the primary variable in focus.) Instead of this automatic determination for each regime, you can choose the same axes to be used throughout the sequence. The automatic choice is recommended in general, as the most dominant variables that are not `fast` or `slow` relative to the primary focused variable of each regime are usually the best to plot on the y-axis. For long animations this may not always be the best situation, however, in order to spot pertinent features in the dynamics.

Regimes with slowly varying variables that act as quasi-static parameters are automatically split up into sub-sequences showing the transitions through the regime, at a resolution of your choice. If you choose to, these sub-sequences can be represented as short AVI movies instead of a large number of on-screen static figures (using movie2avi). A zoom factor fac controls the axes limits of diagrams, relative to the extent of the trajectory segment in each regime. Roughly, the trajectory will occupy 1/fac of the x- and y- axes, although the axes limits will not be extended past the known bounds on the x- and y- axis variables.

9. Domain of Validity / Attractor Estimation

9.1 General guidelines

The type of "attractor estimate" that DSSRT can calculate is currently a form of local `shooting` estimate, with additional epoch consistency checks. It represents a subset of the attractor basin for the reference orbit known to DSSRT. DSSRT will tell you which initial conditions in a variable at a chosen time will lead the system to make qualitatively important changes in the set of Actives within a small tolerance of the epoch times for a known reference orbit. The algorithm is much more efficient than traditional shooting (it only involves one calculation involving all of the data points over the time interval, the rest is based on a binary search), but of course is less accurate. The algorithm is still in development and the behaviour of the algorithm as epochs are traversed is not perfect. An example application of this algorithm is shown in the included PDF slide presentation. The detailed theory behind this algorithm is still being written up.

The `attractor` gives information about the amount of robustness of the orbit in focus to perturbations. Loosely speaking, this is the inverse of the information given by a phase-response curve (PRC) for weakly-coupled oscillators, or a spike-time response curve (STRC) for spiking oscillators. The precise connection between these objects is under investigation. Another important use of the "attractor" is that it represents the domain of validity of an epoch's or regime's reduced model, according to the dominance and time scales set by the epsilon and gamma parameters. In particular, this means that all initial conditions for the system within the computed attractor should be governed accurately by the reduced models.

For the purposes of attractor estimation, it appears to be important not to calculate a functional network using a dominant-scale threshold (epsilon) that is too strict. For neural systems, a good value is between 2.5 and 3.5. The attractor estimation algorithm works best if the time-step in the functional network (set during original integration of orbits, and altered using nOut) is between 0.02 and 0.1 milliseconds. (The time-step is the increment in time observed when stepping through the network in steps of one, using commands 'N' or 'B'.) In this version, inaccurate results can be expected if the time-step is below 1e-5.

There are currently two versions of the AttEst code. One works solely on a transition sequence of epochs. The other works for reduced regimes derived using command `R`. These two versions will be evaluated side by side and possibly the epoch version may be phased out altogether in future releases. The `regimes` version does not do the pseudo-shooting estimation. Instead, for each regime, it estimates the domain of validity of the regime to perturbations in the chosen variable. This is not the same as saying that those perturbations will lead to orbits quantitatively close to the known orbit.

9.2 How to calculate an "attractor estimate"

With an appropriate functional network calculated, focus on the external variable of interest (using command '5'), and set the marks to delimit the time interval of interest. Setting these at epoch transition boundaries for the focused variable gives the best results.

9.2.1 AttEst for epochs only

The next step is to create a transition sequence (command 'C') and agree to focus it on the selected variable. In doing this, it tends to be irrelevant whether the 'ignore order of Acts/Pots' toggle is set or not (command '7'). Red dots will appear above the time bar to show the epoch transition times. (If a transition sequence is loaded from a file then this display of epoch times will not occur.) If a new transition sequence is calculated then the epoch transition times will be refreshed in the diagram. If at any time you need to see which transition sequence is being used to display the epoch transitions, then press 'X' to view the sequences held in memory.

Provided there are more than two epochs calculated in the chosen interval, activate the attractor estimation algorithm using the command 'A'. Select the epoch version of the algorithm. Use the saved values in the parameters files for the supplied demos to get an idea of how what to use for neural systems. Be careful when loading saved attractor estimate parameters from a parameters .par file in case that file specifies a different epsilon threshold value to the one used to create the transition sequence. The correct threshold parameter will be needed by the attractor estimate algorithm.

9.2.2 AttEst for reduced regimes

A set of reduced regimes must first be calculated using command `R` (see Section 8. Reduced Dynamical Regimes) on a fresh transition sequence. Once this has been done follow the above steps in Section 9.2.1 AttEst for epochs only, except select the regime version when prompted (default). In the current version, if the time-step is too small glitches in the "envelope" may occur. I'm looking into why this happens.

Another source of glitches in the results is when your numerical data isn't close enough to the true limit cycle -- if your first regime starts on an epoch other than the first in the transition sequence used. If in doubt, see which epoch RegimeDet wants to start Regime 1 at, and recalculate the transition sequence from that point so that Regime 1 and Epoch 1 coincide.

9.3.3 Interpreting the results

On completion, a plot of the attractor estimate will appear. The original orbit for the focused variable is shown in blue. Vertical red lines mark the ends of epochs that were originally calculated for the functional network. Not shown are the additional epoch boundaries that were used to bolster the temporal resolution of the set of samples used in the attractor estimate calculation. However, they do determine the envelope of the cross-sections, shown using green lines. The green lines provide a linear interpolation of the attractor basin between the sample points.

Multiple attractor estimates (e.g. for different variables) can be held in memory. Any attractor estimate can be redisplayed using the command '&'.

9.3 'Bolstering' the original epoch/regime set

An accurate attractor estimate relies on frequent cross-sections along the orbit at which to analyse the transverse contraction or expansion of perturbations. For many orbits, epoch transitions can be too sparse in some parts, and accuracy is greatly degraded. 'Bolstering' is a simple way to pad out the original epoch/regime set with additional sample points, in order to provide a sufficient resolution of sample cross-sections. These should be added at a higher density over parts of the orbit which are changing most quickly. This is done automatically by AttEst (for both the epoch and regime versions of the algorithm), and is controlled by three simple parameters, providing two resolutions. The parameters that control the bolstering of the original epoch set are:

· The search resolution. This is the time-step granularity used to search the variable's orbit for sudden changes that demand bolster intervals. It also defines the high-resolution granularity of the additional intervals. It must be a positive integer. This integer multiplied by your data's time-step should generally be around 0.05 to 0.1 ms for spiking neural models (Hodgkin-Huxley type spikes).

· The low resolution multiple is the number of search-resolution time units that will constitute a low-resolution granularity of bolster interval. These are the typical-sized intervals added to pad out parts of the orbit that are sparse in real epoch/regime transitions, when the orbit has a small time-derivative. The value must be a positive integer. The low resolution should work out to be around 2ms, for regular spiking neural models. (Slow systems may need both high and low resolutions decreased to improve performance.)

· The derivative threshold measures the relative change in the focused variable necessary to activate high-resolution granularity bolster intervals. High-resolution is activated when the change in the focused variable over one search step, relative to the size of that variable's declared bounds (see the .cfg file) is greater than derivThresh * dt * searchRes, where dt is the absolute time-step in the original variable data. The value must be a positive real number.

· The inner layer resolution is a parameter that cuts off the search algorithm when the variable being perturbed is this close to the original, reference orbit. It plays a similar role to the absolute tolerance in a numerical integrator. Its value should generally be kept very small (and positive) in order to retain maximum accuracy: 0.1 for voltage variables appears to maintain accuracy, and 0.01 for gating variables. The default value is set according to the focused variable's specified bounds automatically, and generally does not need to be changed by the user.

All of the above parameters are user-definable when activating AttEst using the 'A' command, but can be left at the default values in most cases (especially for neural systems). These values may need to be changed, depending on the original sampling time step of the variable data, and the type of ODE system, among other things. At this time, be guided by the default values, and experiment yourself. I hope to write a more detailed account of these parameters for a later version.

In particular, if you see a very irregular profile in an attractor estimate (unexpected sudden changes between consecutive sample points) then you may not have enough bolster intervals: try reducing the search resolutions and possibly the derivative threshold. Correct attractor estimates tend to change smoothly compared to the sample resolution, and tend to change quickly only in accompaniment with obvious large changes in the variables. At this time, for some variables AttEst may not be tunable to always avoid irregularity in parts of the attractor estimate. A fix for this is in development.

In a final point, lowering the low-res multiple parameter too much (for an absolute time-resolution below 1ms) seems to give non-intuitive results. In particular, it makes the attractor estimate substantially decrease in size (as well as greatly increasing the computation time). This may be a 'soft' bug in the algorithm, and I'm looking into it.

9.4 How DSSRT deems what are qualitatively important changes in Acts

For a variable of focus, V, the "fast V-dependent variables" are all variables which (a) have an associated differential equation, and (b) have a time scale for the current epoch/regime that is very fast in comparison to that of V itself. The threshold ratio of the time scales in (b) is given by an internal (non-user definable) parameter. This is set at 2.5 by default, which appears to be ideal for neural models.

AttEst maps a target interval in V back through each epoch/regime included in the time interval of interest, using inverse flow maps generated from the reference orbit and the differential equations of the system (see PDF theory paper for a detailed explanation of this). For all V values in the mapped intervals, their associated Acts set is generated and compared to that of the reference system at that event time. Essentially, AttEst will tolerate discrepancies between the sets when those involve non-V dependent variables or non-fast V-dependent variables, with certain caveats. The exact algorithm will be detailed here in a future release, but can be found in the function CompareActs() within AttEst.m. When V-dependent variables are involved, AttEst decides whether it is more accurate to use the original value of that variable from the reference orbit, under the perturbation in V, or the asymptotic target value if it is deemed "fast".

10. Troubleshooting and other practical tips

· Pressing the Escape key during Functional Network calculation (command 'U'), Attractor Estimates (command 'A'), or most other long calculations, will interrupt and terminate those operations, and return to the main window. This possibility is notified in the console at the beginning of the operation.

· If the KEEP_LOG option (in DSSRT_useroptions.m) is used, periodically delete the session history file DSSRT_History.txt because it will grow unchecked, accumulating with each session.

· Attractor estimate calculations run much slower with the verbose switch on. Unless you are following the detailed behaviour of this function during its calculations, switch it off with the command '?' before running AttEst. Or you can turn verbosity off by default by setting the switch in the user options file DSSRT_useroptions.m.

· For commands that require many dialog boxes to set up parameter values, note that most have highlighted defaults that can be selected just by pressing the space bar. This way you can flick through the default setup (or previous setup, if changes were made) without making many slow mouse clicks each time.

· On occasion you may find that the drawing window for StepNet messes up, possibly in conjunction with another Matlab figure window being misdrawn. This is usually because the StepNet main window is re-selected at a crucial time before a different figure is updated. While the program attempts to avoid this occurrence, it appears to be untrappable in some cases. To try to recover from this, press '0' to redraw the main window. You may have to close the corrupted figures and re-display them (the most recently computed attractor estimates can be redisplayed using '&').

· In my experience with DSSRT running on Windows NT machines, compiling any of the m-files to mex-files except those already listed in DSSRT.m will actually degrade performance, not improve it.

· DSSRT can be fussy when trying to change the data file in the middle of a session, because there may be some residual bugs associated with such operations. If in doubt, restart DSSRT, but please let me know about particular problems you have found so that I can address them in later versions.

11. Future applications, contributions

If you have an ODE system you would like to analyze using the dominant-scale technique and you would like help setting up your network for use with DSSRT then please get in touch with me. I am very interested in working on new collaborations, both within neuroscience and in other areas of applied mathematics.

12. Glossary of new terms

(Coming soon ... for now this is just a lexicon. Use your text editor's search facility to track down the definitions in the main text or in the PDF documentation)

dynamic variable

passive variable (non-dynamic variable)

input (to a differential equation)

input term

simple input (non cross-multiplying variables in input term)

primary variable (for cross-multiplying variables in non-simple inputs)

auxiliary variable (for cross-multiplying variables in non-simple inputs)

observable variable

external variable (externally-observable)

internal variable (not externally observable)

influence scale

dominant scale threshold, epsilon (for inputs)

time scale

time scale threshold, gamma (for variables)

small variation threshold (of a variable)

dominance

dominance strength (of inputs)

asymptotic target value (of a dynamic variable)

focused variable

active (variable)

potential (potentially active variable, potentiated)

inactive (variable)

active set (`Acts`)

potentials set (`Pots`)

functional network

node

link

state map

event

epoch

(event) transition sequence (TS)

fresh transition sequences

regime

slow, fast (characteristic) time scale of a variable's response

near-constant variable