1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104

<html>
<head>
<title>Direct Method with TimeDependent Propensities</title>
</head>
<body>
<h1>Direct Method with TimeDependent Propensities</h1>
<p>
Cain has an implementation of the direct method
for systems of reactions that have timedependent propensities.
To use it, select the "Time Inhomogenous" time dependence
category in the method editor.
Let α(t) be the sum of the propensities. If each of the
propensities was approximately constant on the time scale of 1/α(t),
which is the average time to the next reaction,
then an approximate solution method could treat them as if they were
actually constant. Of course one would need to evaluate all of the
propensities at each step. If any of the propensities varied significantly
on that time scale then we would need to account for this behavior.
In the following exposition we will assume that no propensities become
zero during a step.
</p>
<p>
Consider the exponential distribution with rate parameter λ.
The probability density function is λ e<sup>λ t</sup>;
the mean is 1/λ. Let E be an exponential deviate with unit rate
constant. We can obtain an exponential deviate with rate constant λ
simply by dividing by λ. Now consider the case that the rate
parameter is not constant. A exponential deviate is T where
∫<sub>0</sub><sup>T</sup> λ(t)dt = E. Note that for constant
λ this equation reduces to λ T = E.
</p>
<p>
Recall that when using the direct method one uses exponential deviates
to determine when reactions fire. To determine the time to the next reaction
we generate a unit exponential deviate and then divide that by the sum
of the propensities α. This gives us an exponential deviate with
rate parameter α. Now consider a system of reactions in which the
reaction propensities are functions of time. In order to determine
the time to the next reaction we need to generate a unit exponential
deviate E and then numerically solve
∫<sub>t</sub><sup>t+T</sup> α(x)dx = E for T.
</p>
<p>
To solve for T we can numerically integrate α(t). Below is a simple
algorithm for this.
<pre>
T = 0
while α(t+T) Δt < E:
E = α(t+T) Δt
T += Δt
T += E / α(T)
</pre>
</p>
<p>
You might recognize the above algorithm as the
<a href="http://en.wikipedia.org/wiki/Forward_Euler_method">
forward Euler method</a>, the simplest method for integrating ordinary
differential equations. The accuracy of this method depends on
Δt. There are more accurate methods of numerically integrating
α(t). The
<a href="http://en.wikipedia.org/wiki/Midpoint_method">
midpoint method</a> and the
<a href="http://en.wikipedia.org/wiki/RungeKutta_methods">
fourthorder RungeKutta method</a> are good options.
</p>
<p>
So now we know how to determine when the next reaction fires, but how do we
determine which reaction fires? To do this, we integrate each of the
reaction propensities:
pmf<sub>i</sub> = ∫<sub>t</sub><sup>t+T</sup> a<sub>i</sub>(x)dx. To
select a reaction we draw a discrete deviate with this weighted probability
mass function. Below we use the forward Euler method to calculate the
time step T and the probability mass function pmf used to pick a reaction.
We assume that Δt has been initialized to an appropriate value.
<pre>
s = 0
for i in 1..N:
pmf<sub>i</sub> = 0
p<sub>i</sub> = a<sub>i</sub>(t)
s += p<sub>i</sub>
T = 0
while s Δt < E:
E = s Δt
T += Δt
for i in 1..N:
pmf<sub>i</sub> += p<sub>i</sub> Δt
p<sub>i</sub> = a<sub>i</sub>(t+T)
s += pmf<sub>i</sub>
Δt = E / s
T += Δt
for i in 1..N:
pmf<sub>i</sub> += p<sub>i</sub> Δt
</pre>
</p>
</body>
</html>
