/* MODEL NUMBER:
// MODEL NUMBER: 0329
MODEL NAME: Brusselator
SHORT DESCRIPTION: The Brusselator (combination of Brussels and Oscillator) equations of
I. Prigogine and R. Lefever are solved both continuously and stochastically.
*/
import nsrunit;
unit conversion on;
math Brusselator { // ODE
realDomain t sec; t.min=0; t.max=233; t.delta=0.1;
real A = 5 mM, B= 30.0 mM;
real X0= 30 mM, Y0 = 30 mM;
real X(t) mM, Y(t) mM;
when(t=t.min) {X=X0; Y=Y0;}
real k1 = 0.1 sec^(-1);
real k2 = 0.001 mM^(-2)*sec^(-1);
real k3 = 0.01 mM^(-1)*sec^(-1);
real k4 = 0.1 sec^(-1);
Y:t =-k2*X*X*Y+k3*B*X;
X:t = k1*A+k2*X*X*Y-k3*B*X-k4*X;
}
/* DIAGRAM:
Four chemical reactions are being modeled:
Reaction Rate Constant
A -> X (k1)
2X + Y -> 3X (k2)
B + X -> Y + C (k3)
X -> D (k4)
DETAILED DESCRIPTION: The Brusselator equations form an oscillating system that
approaches a limit cycle. For an example of oscillating systems watch a video
of the Briggs-Rauscher Oscillating Reaction (several videos available on the
internet).
The same equations are solved stochastically using the Gillespie method.
A and B are fixed, X and Y vary, C and D are not calculated.
NOTE: to increase the number of molecules by a factor
of 10, divide k1 by 10, k2 by 1000, k3 by 100, and k4 by 10.
Be sure to make the corresponding change in the ODE model.
Compute lambda1, lambda2, ... lambda4 as shown in the code
for model StochasticBrusselator.
Sum the lambda's to calculate lambda.
Create a random exponential distribution by taking the natural
log of uniform random numbers.
When multiplied by -1/lambda it becomes the time step when
another reaction takes place. The uneven time steps are
accumulated in tau.
A second uniform random number for which reaction takes place
is determined by the probability of the reaction which is the
ratio of its lambda to the sum of all the other lambdas.
We Calculate P1=lambda1/lambda, P2=(lambda1+lambda2)/lambda, ...
For the reaction selected, the appropriate pools are
incremented or decremented in whole units.
KEY WORDS:
Brusselator, stochastic, Gillespie, oscillating reaction, equilibria,
limit cycle
REFERENCES:
I. Prigogine and R. Lefever, Symmetry Breaking Instabilities in Dissipative Systems. II,
J. Chem. Phys. 48, 1695 (1968).
D.T. Gillespie, Exact stochastic simulation of coupled chemical reactions,
J. Phys. Chem., 1977, 81 (25), pp 2340–2361.
REVISION HISTORY:
JSim SOFTWARE COPYRIGHT AND REQUEST FOR ACKNOWLEDGMENT OF USE:
JSim software was developed with support from NIH grants HL088516,
and HL073598. Please cite these grants in any publication for which
this software is used and send one reprint of published abstracts or
articles to the address given below. Academic use is unrestricted.
Software may be copied so long as this copyright notice is included.
Copyright (C) 1999-2011 University of Washington.
Contact Information:
The National Simulation Resource,
Director J. B. Bassingthwaighte,
Department of Bioengineering,
University of Washington, Seattle, WA
98195-5061
*/
import nsrunit;
unit conversion off;
math stochasticBrusselator {
// INDEPENDENT VARIABLE
realDomain t sec; t.min=0; t.max=104; t.delta=0.1;
// PARAMETERS
real k1 = 0.1 sec^(-1);
real k2 = 0.001 sec^(-1);
real k3 = 0.01 sec^(-1);
real k4 = 0.1 sec^(-1);
// THE lambda_i's are functions of the rate constants and concentrations
real lambda1(t) sec^(-1);
real lambda2(t) sec^(-1);
real lambda3(t) sec^(-1);
real lambda4(t) sec^(-1);
// DEPENDENT VARIABLES
realState tau(t) sec; // time when a reaction occurs
when(t=t.min) {tau=0;}
realState X(t), Y(t);
real A=5, B=30, X0=30,Y0=30;
when(t=t.min) {X=X0; Y=Y0;}
// CALCULATE THE lambda_i's
lambda1=A*k1;
lambda2=X*(X-1)*Y*k2;
lambda3=B*X*k3;
lambda4=X*k4;
real lambda(t)=lambda1+lambda2+lambda3+lambda4;
// TWO UNIFORM RANDOM NUMBERS
// The first random number will determine the time of the next reaction
// The second random number determines which reaction
real uni1(t) = random();
real uni2(t) = random();
real expDist(t);
expDist = -1.0/lambda*ln(uni1);
event(t>t.min) tau=tau+expDist;
real P1(t), P2(t), P3(t), P4(t);
P1=lambda1/lambda;
P2=(lambda1+lambda2)/lambda;
P3=(lambda1+lambda2+lambda3)/lambda;
P4=(lambda1+lambda2+lambda3+lambda4)/lambda;
event(uni2<P1) { X=X+1; }
event(uni2>=P1 and uni2<=P2) { X=X+1; Y=Y-1; }
event(uni2>P2 and uni2<=P3) { X=X-1; Y=Y+1; }
event(uni2>P3) { X=X-1; }
}
/*
The Brusselator equations form an oscillating system that
approaches a limit cycle. For an example of oscillating systems watch a video
of the Briggs-Rauscher Oscillating Reaction (several videos available on the
internet).
The same equations are solved stochastically using the Gillespie method.
A and B are fixed, X and Y vary, C and D are not calculated.
NOTE: to increase the number of molecules by a factor
of 10, divide k1 by 10, k2 by 1000, k3 by 100, and k4 by 10.
Be sure to make the corresponding change in the ODE model.
Compute lambda1, lambda2, ... lambda4 as shown in the code
for model StochasticBrusselator.
Sum the lambda's to calculate lambda.
Create a random exponential distribution by taking the natural
log of uniform random numbers.
When multiplied by -1/lambda it becomes the time step when
another reaction takes place. The uneven time steps are
accumulated in tau.
A second uniform random number for which reaction takes place
is determined by the probability of the reaction which is the
ratio of its lambda to the sum of all the other lambdas.
We Calculate P1=lambda1/lambda, P2=(lambda1+lambda2)/lambda, ...
For the reaction selected, the appropriate pools are
incremented or decremented in whole units.
*/