based on the
KINSIM User's Manual
by Bruce A. Barshop
Washington University
Biological Chemistry Computing Facility
Links to external sites may not display from within the program Help window. Use your Web browser to see the original manual at http://bililite.com/tenua/Tenua-manual.html.
Chapter 1: Short Introduction to Using Tenua
Chapter 3: Mathematical Matters
Chapter 4: Writing a Mechanism Description
The Output Expressions Section
Chapter 5: Writing a Data File
Chapter 6: Writing a Variable Values File
Chapter 8: Working with Data in Scripts
Welcome to Tenua, the kinetic simulation program for Java, inspired by my earlier HopKINSIM, which was based on Barshop, Wren and Freiden's KINSIM for the VAX. Tenua simulates chemical reactions by deriving differential equations for the concentrations of chemical species from descriptions of the reactions, then numerically solving those differential equations.
This manual assumes you know what it means to get a differential equation from a reaction. If not, then there are some good tutorials available on the web from various sources, such as the University of Delaware, Purdue University, and chemistrycoach.com.
The name Tenua is from the Hebrew word for "movement," which is the closest I could get to "kinetics."
There are many other kinetic simulation programs available. The grandaddy of them all is KINSIM, which also includes a fitting program to fit mechanism parameters to real data. Unfortunately, it is a DOS program, and has not been updated since 1997. The original program was described in BA Barshop, RF Wren, and C Frieden (1983). Analysis of Numerical Methods for Computer Simulation of Kinetic Processes: Development of KINSIM-A flexible, portable system. Analytical Biochem. 130, 134-145. A more modern version is available from KinTec, under the name KinTecSim, that is designed to work with the data generated by their stopped-flow machines. An excellent numerical-integration kinetics program is Dynafit, from BioKin. It allows for detailed control of the mechanism and all the numerical parameters, with very flexible curve fitting and allows for perturbation-from-equilibrium experiments, which would be hard to do with any of the KINSIM-based programs. Many of the web-based tutorials also include simple kinetic simulation programs.
All these use numerical integration of derived differential equations to solve the simulation. A different approach is to directly simulate the collisions and reactions of individual molecules or ensembles of molecules. This can often solve problems not well-described by differential equations, and allow for variations in things like reaction volume and pressure. A good visualization of a bimolecular reaction at the molecular level is here. A very powerful and flexible program from IBM that uses this approach is CKS.
Tenua is "e-postcardware;" it is freely distributable (see the license) as long as all these files are included with it, but I would appreciate an emailed picture of your (hopefully scenic) environs to hang on my wall (and to give me an address to send upgrades to). Please send your photos and bug reports (and lists of desired new features and any cash you can spare) to:
Danny Wachsstock
d.wachss@prodigy.net
The program and the source code released are subject to the Sun Public License Version 1.0 (the License); you may not use this program except in compliance with the License. A copy of the License is available from this site or from sun.
Tenua is a Java 1.4 program, so it should run on Windows, Macintosh and Linux systems without modification. It is distributed as a Java archive (.jar) file that can be run by just double-clicking on the icon. For Macintosh users, the program may look unfamiliar, with the menu bar at the top of each window instead of the top of the screen.
Using Tenua consists of 3 steps:
A + B <-> C // the double slash indicates a comment
// This is a bimolecular reaction, with the double headed arrow
// sign representing a kinetic reaction
// (it's a less-than sign, dash, greater-than sign)
startTime), the ending time (endTime),
and the intervals at which to produce results (timeStep). There is also a variable called
time, but the program will set that, starting at startTime
and successively adding timeStep until it exceeds endTime.
You can also set epsilon, which is the relative accuracy of the numerical solution.
The smaller epsilon is, the more accurate the simulation will be,
but it will be that much slower. A reasonable number is 1e-4, or 0.01%.k(+1)
being the forward rate
constant for the first reaction and k(-2) being the backward rate constant for the second reaction.
To the left of that is the mechanism parameters, which are any variables that are not one of the above.
a=1; b=2; k(+1)=3; endTime=4;. Spacing and separating
lines will not affect it, but there must be a semicolon after every statement.
See Chapter 6 for more details
The menu items are mostly self-explanatory. The toolbar icons are all the same as the corresponding menu items, but are there for easy access.
var1 = 1; var2 = 2; k(+1) = 1e-6; etc.
See Chapter 6 for details.
If there are errors, an error message window will be shown.
1;2;3;4 adds 4 new times to the data if those
times did not already exist.
Select the Table tab and double-click a given cell to enter a value.
First;Second;Third adds 3 new columns to the data, if those names did
not already exist.
The text boxes at the left of the graph control the scale of the axes. If the Auto-adjust Axes is checked, then the scale will be expanded as data are added to show all of the data. If you want to zoom in on a section of the graph, uncheck this box, then enter the desired scale in the text boxes. Press the Optimize Axes to shrink or expand the axes to show all of the data, such that the plots fill the graph.
The Reset button sets all the variables to their default values.
This tab is simply a text editor where you can record notes and save or print them. You can add notes from within a script using the note commands.
Understanding how Tenua works is key to getting it to work well. The compiler takes a mechanism description like
A <-> B <-> C
And turns it into a series of calculations like
rate(+1)=k(+1)*A;
rate(-1)=k(-1)*B;
rate(+2)=k(+2)*B;
rate(-2)=k(-2)*C;
rate(A) = -rate(+1)+rate(-1);
rate(B) = rate(+1)-rate(-1)-rate(+2)+rate(-2);
rate(C) = rate(+2)-rate(-2);
It then solves the ordinary differential equations represented by the
rate(A), rate(B), and rate(C) equations. The
algorithms for solving the equations are taken from the book,
Numerical Recipes in C++, chapter 16.
The initial values are generally taken from the Initial Variable Values tab, but can be overridden in the script section of the mechanism description (see Chapter 4). The Initial Variable Values tab gets these values from initialization statements in the mechanism, like:
A: 1; k(+1): 10;
If there is no initialization statement, then the program uses the most recent value for those variables in some other mechanism, or zero otherwise.
Note that these calculations are all unitless. The program assumes everything
is in the same units. If the concentrations represent mM and startTime
and endTime are in s, then the rates are mM/s.
The closer the numbers are all to 1, the more efficient the differential equation solver. So choose your units appropriately! If the rate constants differ by more than a few orders of magnitude, then there may be no way to get all the rates and concentrations to be on the order on 1. If this is a problem, see Chapter 7: Stiff Equations.
The compiler also puts all the other calculations (as in the output section) into a list of calculations to be evaluated at each time step. I'll call these the output calculations.
The simulation runs by setting time=startTime,
evaluating the output calculations and inserting the values into the table.
It then subdivides timeStep into a number of substeps and for each substep,
evaluates all the rate(x)s for each chemical x
and calculates each x as
x = x + rate(x)*substep.
It then sets time = time + timeStep,
evaluates the output calculations and inserts the values into the table.
It continues this cycle until time is greater than
or equal to endTime.
The number of substeps and how to calculate them depends on the
algorithm used for the differential equation solver and the value of epsilon.
If timeStep is 0, then the program tries to select
values for timestep that will give the right accuracy for the
given epsilon, which may be all you need.
If the graph looks jerkier than you want, then you will need to set timeStep.
The textual mechanism description consists of three major sections: the first section contains the chemical equations, initialization statements and intermediate calculations; the second contains the output equations; the third contains the "script". The chemical equations define the reaction scheme and the output equations define the variables to be graphed during the simulation. The script describes exactly how to run the simulation.
The syntax is based closely on Java's syntax, so every line or statement must end with a semicolon. The formal grammar is in a separate file; this is less rigorous discussion.
Comments: A double slash (//) in the text file indicates that
the remainder of that line is a comment to be ignored.
A comment in the middle of a line is enclosed in
/* ... */.
The first section of a mechanism description is the reactions section. It contains chemical equations, initialization statements and intermediate calculations, each ended with a semicolon. They may be freely intermixed.
These equations are entered in standard chemical format.
Each species in the mechanism is represented by a unique name that
starts with a letter and has any combination of letters, numbers
and the symbols ', $ or _.
Numbers preceding a species name signify the stoichiometry in the reaction.
For example:
Note that numbers within names have no meaning per se. ThusAmay represent a monomeric species
A'may represent an activated form ofA
AAorA2may represent the dimer of that species
2A2represents two molecules ofA2
4AB3O4represents four molecules namedAB3O4
A2 is not necessarily a dimer of A,
and the program will not enforce any mass conservation based on the
names of species. The only mass conservation is based on the chemical
reactions as written.
Chemical steps in the mechanism are represented by double-headed
arrows (<->). This signifies a step which is rate-governed,
that is, a step with both a forward and reverse rate constant.
Consecutive steps may be concatenated:
E + S <-> ES <-> EP <-> E + P;
Non-consecutive steps may be entered on one line, separated by semicolons:
E + S <-> ES ; E + I <-> EI ;
Any variable can be given an initial value in the mechanism description.
These include chemical species (names that appear in the chemical equations),
rate constants (which have the form k(+1) for the forward rate
constant for the first reaction and k(-2) for the backward rate constant
for the second reaction), time constants (things like endTime
and timeStep) and parameters (names that do not appear anywhere
else, that are used for intermediate calculations).
An initialization statement uses a colon to indicate the value; it has the form:
name : 123.45 ;
Statements can be concatenated, so the following is legal:
k(+1) : k(-1) : k(+2) : 1;
Rate constants can use the shorthand of k(+) for the forward rate constant of
the most recently defined reaction, and k(-) for the backward rate constant.
Statements can include arbitrarily complex expressions and use previously defined variables, so the following are legal:
A : 1+2^4;
B : A * 2;
See the next section for more details on expressions.
Initialization statements are evaluated when the mechanism description is compiled, not each time the mechanism is run. This means that you can change the initial value of a variable in the Initial Variable Values tab even if it was initialized in the mechanism description.
There is a special form of the initialization statement to define constants:
final pi : 3.14159 ;
defines a constant called pi that cannot be changed in the
Initial Variable Values tab. Some constants are predefined; they are
final avogadro : 6.0221367e23;
final boltzmann : 1.380659e-23; // in J/K
final gasConstant : 8.31451; // in J/(mol*K)
final infinity: Double.POSITIVE_INFINITY; // represents an infinitely large number; this can be useful in curve fitting
pi is not one of the predefined constants; it was used as an example of defining
your own constants in a mechanism.
Calculation expressions use an equal sign to set a value at each time step during the calculation. Thus:
name = 1 ;forces
name to be 1 at each time step. This is different from intialization, where
name : 1 ;sets
name to be 1 at the beginning of the simulation,
and lets you change it before simulating using the Initial Variable Values tab.
The expressions are specified in a Java-style algebraic notation, with the following considerations:
Names of species represent their concentrations at the current
time in the simulation. Names not previously used refer to parameters.
The variables time, startTime, endTime, timeStep
and epsilon are set by the program and should not be assigned to (but they should be
initialized).
Numbers are IEEE double-precision numbers. To use scientific notation, use e,
as is 6.02e23.
The mathematical operators defined are:
+ - * / % x % y is the remainder of x/y,
so 5 % 2 is 1.
^ 2 ^ 3 is 8. This is not a Java
operator, so it can be only used in the reaction and output expressions
sections, not in the script.
1; if false, it evaluates to 0.
These are defined in Java, but differently from the way they are used here. In the script section,
they have their Java definitions.
== 2==2 is 1;
2==3 is 0. The double-equals sign is from Java; it
lets the compiler distinguish between assignment and the equals operator.
!= < <= > >= & x & y is 1
if x and y are both non-zero,
and is 0 otherwise.
| x | y is 0
if x and y are both zero,
and is 1 otherwise.
? x ? y is equals to the value of y
if x is non-zero, and is equal to 0 otherwise.
Mathematical operators have the precedence:
^higher than*,/,%higher than+,-higher than==,!=,<,<=,>,>=higher than&,|higher than?
Parentheses can be used to resolve ambiguous expressions. For example:
A+B*C adds A to the product of B and C and
(A+ -B)*C multiplies C
by the sum of A and negative (B).
Four functions are defined: -x is negative x, like you would expect;
!x represents the Boolean not operator (if x is non-zero, then
!x is 0; if x is zero, then !x
is 1);
ln(x) is the natural logarithm of x; and exp(x)
is for Euler's number e^x.
For example, k(+1) = A*exp(-Ea/(boltzmann*temp)); lets you set the rate constant
of the first reaction using the Arrhenius equation, where A,
Ea and temp are parameters and boltzmann is a predefined
constant.
Functions that start with @ (the at sign) represent Java functions
that are defined outside the mechanism. You can use all of the
java.lang.Math library
or any other java library. For example,
@Math.sin(x) returns the sine of x,
and @System.currentTimeMillis() returns the real time in milliseconds
from midnight, January 1, 1970 (this can let you time how long a simulation takes).
You have to make sure that the function you are calling returns a number;
the simulation will show an error dialog if it does not. The error message will be
something like "Index out of bounds."
You can define your own functions in the script section;
see the script section. To call them,
use @script.functionName(x).
You can also explicitly change the way the rates of change of the chemical reactions are set up. To recap Chapter 3, a chemical equation of:
a <-> b;becomes:
rate(+1) = k(+1)*a;
rate(-1) = k(-1)*b;
rate(a) = -rate(+1)+rate(-1);
rate(b) = rate(+1)-rate(-1);
If this is not right, then you can explicitly assign rate(+n)
and rate(x) for any reaction n or chemical x.
For example, if a is a source or a sink
(i.e., its concentration should not change), then write:
a <-> b;
rate(a) = 0;
or if the forward rate is some more complex function of a, then write:
a <-> b;
rate(+) = @script.function(a);
Note that rate(+) is a shorthand for the forward rate for
the most recently defined reaction.
Note also that using a = some_constant; would not define
a constant-concentration source, because simple intermediate expressions
are only calculated at the end of each time step, and the source needs
to be constant during the differential equation calculations as well.
Similarly, time is only updated at the end of each time step.
So if a rate depends on time, you need to define a new parameter (say,
call it t), and use:
t: startTime; rate(t) = 1;then use
t in your expression for that rate.
See the Numerical Recipes book,
equation 16.6.19 and the discussion there.
This means that you can use Tenua for solving arbitrary ordinary differential equations. To solve dx/dt=sin(x^2) , use:
rate(x)=@sin(x^2);
x : 0; // initial conditions
The program will assume that x here is a chemical that was
not in any reaction yet.
After the reactions section, you can have a section that defines what
expressions are to be output to the table and graph. This section
starts with *output, exactly like that, then a
list of semicolon-terminated expressions (see the expression section
for a description of expressions). The name of the output expression
is the expression itself, unless you assign it to a variable.
In that case, the name of the output is the name of the variable. Thus:
*output
a; // output the value of a, labeled "a"
sample=k(+1)*b; // output the value of k(+1)*b, labeled "sample"
If there is no output expression section, then all the chemical species concentrations will be output. So a mechanism description of:
A<->B<->C;is the same as:
A <-> B <-> C;
*output A; B; C;
The script is an optional section that starts with *script,
exactly like that, then a series of commands that each ends with a
semicolon. The possible commands are:
go;
The go command is the most important one; it actually
runs the simulation once with the currently-defined initial variable
values. If the *script section is not present, the
default script is simply
*script
go;
If you want to run the simulation twice with the same variables
(I don't know why you would), use
*script
go;
go;
rename(n,newName);
The rename command changes the name of the nth
output to newName for next time go is executed.
This is important because every time go is executed, the
new outputs overwrite any old data column with the same name. So, to
run a simulation twice but keep the data for both runs, use
a <-> b;
*output
a; b;
*script
rename (1, "a-first run");
rename (2, "b-first run");
go;
rename (1, "a-second run");
rename (2, "b-second run");
go;
The new name in the rename command must be a string
surrounded by quotes or a number, but the compiler is smart enough
to combine strings and numbers (using +). So
rename (1, "k(+1) = " + k(+1)) where k(+1)
is 1.5, would set the name of the first output to "k(+1) = 1.5"
Displays the string in a dialog box and gives you the option of stopping the simulation or continuing. The name comes from Visual Basic. A part of an example:
go;
msgbox("Do you want to do it again?");
go;
If the user clicks "Stop" the script is stopped with an error message.
note commands
note(string) adds string to the Notes tab.
noteTime() adds the current date and time to the Notes tab.
note variables var1,var2… adds a series of lines
to the Notes tab, with the names and values of the named mechanism variables.
This can be useful for automated curve fitting to ouput the fitted
parameters.Example:
A<->B;
k(+1) : 1;
A : 10;
*script
note("Starting simulation");
noteTime();
go;
note("Ended simulation");
noteTime();
note variables A, k(+1);
produces the following output on the Notes tab (for the time I ran it):
Starting simulation
// Aug 12, 2005 1:39:39 AM
Ended simulation
// Aug 12, 2005 1:39:39 AM
A=10.0;
k(+1)=1.0;
The format of the noteTime() and note variables output is
designed to allow you to save the file and use it as a variable values
file.
list(item1,item2,item3...)
i-th item
can be accessed as array-name[i], but the numbering starts
from zero! This is a feature of most modern programming languages.
a = list (10, 20, 30, 40);
msgbox (a[0]); // prints "10"
msgbox (a[3]); // prints "40"
msgbox (a[4]); // generates an error
This is most useful in the Java 1.5 for-each loop:
for (x : list (0.1, 0.2, 0.3, 0.4, 0.5)){
k(+1) = x;
go;
}
This will run the mechanism 5 times, with k(+1) set to each successive
value. The sophisticated programmer might want to write
for (k(+1) : list (0.1, 0.2, 0.3, 0.4, 0.5)){ // DON'T DO THIS
go;
}
but I can't get the for-each loop to work with mechanism variables. You have to use
a variable defined in the script and assign the mechanism variable to that.
a=1; where a is
a variable defined in the mechanism will be interpreted as
changing the initial value of that variable. Note that this overrides
any initialization that was done in the mechanism with
a: 1; or done by setting the value in the
Initial Variable Values tab. So a useful mechanism that
runs more than once might be:
a <-> b;
*output a;
*script
k(+1)=1;
rename (1, "first");
go;
k(+1)=2;
rename (1, "second");
go;
k(+1)=3;
rename (1, "third");
go;
or, using rename in a more sophisticated way,
a <-> b;
*output a;
*script
k(+1)=1;
rename (1, "k(+1) = " + k(+1));
go;
k(+1)=2;
rename (1, "k(+1) = " + k(+1));
go;
k(+1)=3;
rename (1, "k(+1) = " + k(+1));
go;
and this will label the output of each run with the
appropriate value of the rate constant.
Any Java statement can be used in the script. The program uses the BeanShell interpreter to execute the script. This gives you tremendous flexibility, but I cannot teach you all of Java in this manual. For example, this is a simpler version of the above scripts:
a <-> b;
*output a;
*script
for ( k(+1) = 1; k(+1) <= 3; k(+1) += 1){
rename (1, "k(+1) = " + k(+1));
go;
}
For more sophisticated programming, two objects are defined
that give access to the inner workings of the program. You need to download and read the
Javadoc for the program
to fully use them. mechanism is an object that
represents the currently running
mechanism; see the Javadoc for tenua.simulation.Mechanism for the sorts of things
you can do with that. See Chapter 7: Stiff Equations
for one common use of mechanism. dataTable is an object that represents
the table of data; see the the Javadoc for tenua.gui.DataModel for the sorts of things
you can do with that. See Chapter 8: Working with Data in Scripts
for one common use of dataTable.
One sort of Java statement is a method declaration. You can define your own methods and call
them from inside the script, as though the entire script was inside a class declaration. In fact,
this is exactly what the compiler does. If you want to use the method in a calculation in the
mechanism, then the method must return a number, and all the arguments must be numbers. Call it
with @script.methodName(arguments) in the mechanism. See the
sample mechanism descriptions for examples.
The mechanism below represents a simple unimolecular reaction, with all the variables explicitly inititalized.
a <-> b; // the reaction
k(+): 2; // initialization statements
k(-): 1;
a: 10;
b: 0;
startTime: 0;
endTime: 10;
epsilon: 1e-6;
timeStep: 0;
*output // output expressions
a; b;
The mechanism depicted below represents the simple enzyme
reaction catalyzed by fumarase. In this mechanism, E
represents the enzyme, F its substrate fumarate and M the product malate.
// Fumarase mechanism
E + F <-> EF <-> EM <-> E + M; // Reactions
*output // Output expressions follow
EF*X1 + F*X2 - X3; // Two output
EM*X4 + M*X5 - X6; // ...expressions
The output expressions include the names of species in the reaction
scheme as well as output constants. The first output expression will
equal the concentration of the enzyme-substrate complex (EF)
times its extinction coefficient (X1)
plus the concentration of the substrate (F)
times its extinction coefficient (X2) minus some
offset factor (X3). At run-time the values of the factors
may be set for any particular display. For example, if the values
of X1, X2 and X3 are set to
0, 1, and 0 respectively,
the first output will simply equal the concentration of F.
The mechanisms below depict two schemes for bi-bi enzyme reactions. The first one is a compulsory order for substrate addition and the second one has a random addition.
// Ordered Bi-Bi
E + A <-> EA;
EA + B <-> EAB <-> EQ + P;
EQ <-> E + Q;
*output
F1*P;
F2*Q;
F3*E;
F4*A;
F5*B;
F6*EA+F7*EAB;
F8*EQ;
// Random Bi-Bi
E + A <-> EA ; E + B <-> EB;
EA + B = EAB <-> EPQ;
EQ + P <-> EPQ;
E + Q <-> EQ ; E + P <-> EP;
*output
F1*P;
F2*Q;
F3*E;
This example is from the Colby College department of chemistry. It demonstrates how a chemical reaction can sustain oscillations. This requires reactions far from equilibrium that have an autocatalytic step and two stable states.
a+m <-> 2a; k(+): 0.1; k(-): 0;
a+b <-> 2b; k(+): 0.1; k(-): 0;
b <-> p; k(+): 0.05; k(-): 0;
rate(m) = rate(p) = 0; // sink and source
a:1; b:1; m:1; // initialization
epsilon: 1e-4;
startTime: 0;
endTime: 300;
This shows how to use script-defined functions and rate equations
// simulates a solution/dissolution reaction. The disociation constant K = [a][b] is
// k(-1)/k(+1)
a + b <-> ab;
timeStep : .1;
rate(-) = @script.rateFunc (ab);
*output
a; b;
// tenua does not automatically make sure that concentrations are greater than
// zero. With the usual differential equations this is not a problem, since the rate is
// proportional to the concentration, so as the concentration gets smaller, the rate gets
// lower and never forces the concentration below zero. With this zero-order reaction,
// the concentration of ab may go below zero, but only minimally. The smaller timeStep
// or epsilon is, the less of a problem this is. To avoid seeing this problem,
// we output ab only if it is
// greater than zero. The values of a and b will still oscillate slightly around their
// equilibrium values
product = ab > 0 ? ab ;
*script
// the rate of the reverse reaction is zero-order (constant) if there is any ab
// to dissolve, and is zero otherwise
rateFunc (x) { if (x > 0) return k(-1); else return 0; }
go;
Tenua can be used just to solve differential equations, without using chemical reactions.
Just write the equations out with the rate(x)= formalism. For example:
// Solves dx/dt = x^-2 with initial condition x = 1
rate(x) = x^-2;
x : 1;
*output
x;
(3*time+1)^(1/3); // this is the analytic solution, which should be identical
This mechanism shows how you can use System.currentTimeMillis
to see how much real time a simulation takes.
a <-> b;
*script
timeIt(){
timeBegun = System.currentTimeMillis();
go;
timeEnded = System.currentTimeMillis();
msgbox ("Total time = " + (timeEnded-timeBegun));
}
timeIt();
A data file is simply a tab-delimited text file:
Time A B AB 0 .1 .1 0 1 .05 .05 .05 2 .03 .03 .07 3 .02 .02 .08Where the spaces between the columns are all the "tab" character. The first line is the column titles, which are the names used to plot the data. The first column is always interpreted as the time, and the name is ignored.
To generate a tab-delimited text file in Microsoft Excel, create the table in the form above, then select Save As... from the File menu and in the Save as Type drop down box, select Text (tab-delimited).
The variable values file is really just a script in its own file. Anything that can go in a script can go here, so if you want to define a method that can be loaded into any mechanism you can include it in a variable values file and use Load... from the Initial Variable Values tab.
It is designed, however, to load initial variable values, so to initialize the values of the time
constants and the chemicals a and b, create a text file (in Notepad or TeachText or
your favorite text editor) like:
startTime = 0; endTime = 100; timeStep = 1; epsilon = 1e-5; a = 1; b = 0;
If a system of differential equations has a solution with exponential decay rates that differ by many orders of magnitude, then the usual methods of numerically solving it may be unstable. This means that even if one component of the solution is infinitesimally small, it will make the solution as a whole wildly inaccurate unless the step size also made infinitesimally small. For tenua, this means either that the program will fail to converge on a solution, or that it will be unbearably slow. This situation can occur if the rate constants in a mechanism are very different.
The mathematicians call this a "stiff" system of differential equations.
The way to solve this is to change the algorithm that tenua uses to solve the differential equations.
There is a variable called mechanism.solver that is the name of the algorithm. The usual
algorithm is named "normal"; it is the Bulirsch-Stoer algorithm from the Numerical
Recipes book, chapter 16.4. The stiff algorithm is called "stiff"; it is the Bader-Deuflhard
algorithm from chapter 16.6.
To change the algorithm, you have to use a script. Include the line
mechanism.solver="stiff"; before go;. Note that the stiff algorithm will work
for non-stiff problems, it is just slower under those circumstances.
An example of a stiff problem, from Numerical Recipes:
// Numerical Recipes equation 16.6.27-28
rate(y1) = -.013*y1 - 1000*y1*y3;
rate(y2) = -2500*y2*y3;
rate(y3) = -.013*y1 - 1000*y1*y3 - 2500*y2*y3;
y1:1; y2:1; y3:0;
epsilon:1e-3;
endTime:50;
*script
mechanism.solver = "normal";
go;
msgbox ("Finished first run");
mechanism.solver = "stiff";
go;
msgbox ("Finished second run");
When you run this, you will note that the first run is painfully slow. This would be an interesting
mechanism to use the timeIt method from the mechanism examples.
There are several other algorithms implemented. The only other one I think might be useful is the
adaptive-stepsize fourth-order Runge-Kutta method (see Numerical Recipes chapter 16.2) for very
discontinuous functions; use mechanism.solver="Embedded Runge Kutta".
You can refer to the table of data in your scripts. This is necessary for doing automatic curve fitting,
as described in Chapter 9: Curve Fitting, but can be useful in other circumstances.
As a general rule, columns in the data table are always in quotation marks, while mechanism variables
never are. Thus a refers to a variable (chemical or parameter) named a, while
"a" refers to the column named a. Just about any character is legal in a
column name; only the quotation mark itself and the semicolon are not allowed. Useful functions include:
dataTable.addColumn("name of column");
data("name of column").insert(time, value);
data("name of column").firstTime()
Double.POSITIVE_INFINITY if
there are no data points. Double.POSITIVE_INFINITY is a predefined
Java constant representing infinity.
data("name of column").lastTime()
Double.POSITIVE_INFINITY if
there are no data points.
data("name of column").nextTime(t)
for (t = data("col").firstTime(); t <= data("col").firstTime(); t = data("col").nextTime(t)) {...}
data("name of column").value(time)
data("x").value(1) is 0 and data("x").value(2)
is 10, and data("x").value(1.5) is not a data point on the table,
data("x").value(1.5) will return 5.
"col"(t) is the same as
data("col").value(t). "col".firstTime() is the same as
data("col").firstTime(), and similarly for "col".lastTime(),
"col".nextTime(t) and "col".insert(t, value).
for (t = "col".firstTime(); t <= "col".lastTime(); t = "col".nextTime(t)){
msgbox ("col at time=" + t + " is " + "col"(t));
}
Or to add up the sum of squares
of the differences between column 1 and column 2, use:
sum = 0;
for (t = "column 1".firstTime(); t <= "column 1".lastTime(); t = "column 1".nextTime(t)){
diff = "column 1"(t) - "column 2"(t);
sum += diff*diff;
}
go varying command.
a <-> b; k(+): 1; k(-):1;
a: 10;
*output b;
*script
fit ("b", "data");
go varying k(+1), k(-1);
After running this mechanism, the Initial Variable Values tab will display the calculated values for
k(+1) and k(-1).
fit ("simulated data", "real data");
go varying
statement will
find the nonlinear least squares fit.
residual (time, data1, data2)
{ return (data1-data2)*(data1-data2);}. You can define any other residual function
in the script section; for absolute values (Numerical Recipes equation 15.7.8) use
residual (time, data1, data2)
{ return Math.abs(data1-data2);}. For a Lorentzian distribution (Numerical Recipes equation
15.7.10) use residual (time, data1, data2)
{ diff = (data1-data2); return Math.log(1 + diff*diff/2);}. The time is included as
an argument to this function to allow time-dependent variations in the estimated standard deviation
of the errors (the sigma in Numerical Recipes equation 15.1.5)
go varying statementgo varying statement actually runs the algorithm to do the fitting. It consists
of the words go varying followed by a comma-separated list of the mechanism variables
to be adjusted to make the data fit. These variables should have their initial values set to your
best guess as to the right answer. The initial guess sets the scale of the problem and the accuracy of
the result; if you set the initial value for a parameter to 100 and the right answer is 0.001, the
program will likely never get there. You can list any number of the mechanism variables, and
you can use any of them: chemical concentrations, rate constants or parameters. They need to have been
used somehow in the reaction section.list
function described in The Script Section. The data shown were generated with
k(+1) equal to 4.5 with some added random errors. This mechanism starts with a guess of
k(+1) equal to 1.0 and ends up with the solution of k(+1) equal to 4.661.
a<->b;
a: 10;
k(+):1;
k(-):1;
startTime:0;
endTime:10;
timeStep:1;
*output b;
*script
dataColumn = list(0.9740033292751074, 8.046241766507238, 8.236222858487372, 8.442838611406437,
7.806820868958864, 8.07694857577928, 8.302839698933445, 8.380954912747303,
7.802980764099486, 8.48219893144326, 8.745798377531942);
dataTable.addColumn ("data");
for (i = 0; i < dataColumn.length; ++i) {"data".insert (i, dataColumn[i]);}
fit ("b", "data");
go varying k(+1);
msgbox (k(+1));
These show up in the bottom panel of the Editor tab as the compiler checks what you are typing. If you are not done, ignore the errors; the compiler runs once a second.
The errors are displayed with their computer-jargon names, which is something that ought to be fixed in a later version. Some terms may come up:
<EOF> stands for "end of file." It means the mechanism description
ended when the compiler expected something more.
";" and "\n"are internally used to end statements; if they come up, it
means that a statement ended when the compiler expected something more.
A syntax error; something like a + + b that the compiler cannot
interpret. It tries to give you a list of what it thought you meant.
A completely unrecognized character.
A constant defined with final was used in a reaction or an
expression that would change its value.
An attempt to assign a value to an expression, like a+1 = 12.
A syntax error in the script section.
A completely unexpected character was found in the script section.
These come up in a separate window if the program encounters an error condition while executing. The top line is the simple description of the error; the next line is a detailed computer-jargon-laden description of the error, and the bottom section is a list of the program methods that were called when the error was produced; useful for debugging but not to the average user.
This is the most common simulation error. If the program cannot attain the desired
accuracy (as represented by epsilon, then it displays this error. The
solution may be to increase epsilon, reduce timeStep, or
change mechanism.solver. Some problems,
especially with very high stoichiometries, may be unsolvable with this program;
try to change the mechanism to one with more steps but lower stoichiometry, or use
a program like
CKS.
Some other unspecified problem with the simulation. Please contact me, as this probably represents a bug in the program.
Part of the program could not be executed. Please contact me, as this probably represents a bug in the program.
An error that should have been caught by the compiler, like assigning to a constant.
The file or pasted string cannot be interpreted as a list of script statements.
The operating system could not complete the desired operation.
The interpreter found some error in the script. You will need to look at the detail
message. command not found or attempt to resolve method
indicates a method was called that was not defined. no such property
generally indicates that the script tried to do something illegal with a mechanism
variable (the compiler turns all mechanism variables into properties;
time becomes variable("time").value which
BeanShell turns into
mechanism.get("time") or mechanism.set("time",x)).
You will get this error for using some Java constructs like ++ or
for (x:Collection)
with a mechanism variable; something like
for (k(+1) = 1; k(+1) <= 10; k(+1)++) or
for (k(+1): list(1,2,3,4) {...}. Replace those with
for (k(+1) = 1; k(+1) <= 10; k(+1)+=1) or
for (x: list(1,2,3,4) {k(+1)=x;...}.
The last action cannot be undone or redone but Undo or Redo was selected.
© 2007 Daniel Wachsstock, MD