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/Tenuamanual.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 KINSIMA flexible, portable system. Analytical Biochem. 130, 134145. A more modern version is available from KinTec, under the name KinTecSim, that is designed to work with the data generated by their stoppedflow machines. An excellent numericalintegration 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 perturbationfromequilibrium experiments, which would be hard to do with any of the KINSIMbased programs. Many of the webbased 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 welldescribed 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 "epostcardware;" 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 doubleclicking 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 lessthan sign, dash, greaterthan 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 1e4, 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 selfexplanatory. The toolbar icons are all the same as the corresponding menu items, but are there for easy access.
var1 = 1; var2 = 2; k(+1) = 1e6; 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.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 Autoadjust 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 <> CAnd 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. ThusA
may represent a monomeric species
A'
may represent an activated form ofA
AA
orA2
may represent the dimer of that species
2A2
represents two molecules ofA2
4AB3O4
represents 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 doubleheaded
arrows (<>
). This signifies a step which is rategoverned,
that is, a step with both a forward and reverse rate constant.
Consecutive steps may be concatenated:
E + S <> ES <> EP <> E + P;Nonconsecutive 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.380659e23; // 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 Javastyle 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 doubleprecision 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 doubleequals 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 nonzero,
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 nonzero, 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 nonzero, 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 constantconcentration 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 conditionsThe 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 semicolonterminated 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 currentlydefined 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, "afirst run"); rename (2, "bfirst run"); go; rename (1, "asecond run"); rename (2, "bsecond 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 arrayname[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 errorThis is most useful in the Java 1.5 foreach 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 foreach 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: 1e6; 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; // ...expressionsThe 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 enzymesubstrate 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 runtime 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 bibi enzyme reactions. The first one is a compulsory order for substrate addition and the second one has a random addition.
// Ordered BiBi 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 BiBi 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: 1e4; startTime: 0; endTime: 300;
This shows how to use scriptdefined 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 zeroorder 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 zeroorder (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 = " + (timeEndedtimeBegun)); } timeIt();
A data file is simply a tabdelimited 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 tabdelimited 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 (tabdelimited).
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 = 1e5; 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 BulirschStoer algorithm from the Numerical
Recipes book, chapter 16.4. The stiff algorithm is called "stiff"; it is the BaderDeuflhard
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 nonstiff problems, it is just slower under those circumstances.
An example of a stiff problem, from Numerical Recipes:
// Numerical Recipes equation 16.6.2728 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:1e3; 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
adaptivestepsize fourthorder RungeKutta 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 (data1data2)*(data1data2);}
. 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(data1data2);}
. For a Lorentzian distribution (Numerical Recipes equation
15.7.10) use residual (time, data1, data2)
{ diff = (data1data2); return Math.log(1 + diff*diff/2);}
. The time is included as
an argument to this function to allow timedependent 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 commaseparated 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 computerjargon 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 computerjargonladen 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