Cain has an implementation of the direct method for systems of reactions that have time-dependent 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.

Consider the exponential distribution with rate parameter λ.
The probability density function is λ e^{-λ t};
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
∫_{0}^{T} λ(t)dt = E. Note that for constant
λ this equation reduces to λ T = E.

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
∫_{t}^{t+T} α(x)dx = E for T.

To solve for T we can numerically integrate α(t). Below is a simple algorithm for this.

T = 0 while α(t+T) Δt < E: E -= α(t+T) Δt T += Δt T += E / α(T)

You might recognize the above algorithm as the forward Euler method, 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 midpoint method and the fourth-order Runge-Kutta method are good options.

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_{i} = ∫_{t}^{t+T} a_{i}(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.

s = 0 for i in 1..N: pmf_{i}= 0 p_{i}= a_{i}(t) s += p_{i}T = 0 while s Δt < E: E -= s Δt T += Δt for i in 1..N: pmf_{i}+= p_{i}Δt p_{i}= a_{i}(t+T) s += pmf_{i}Δt = E / s T += Δt for i in 1..N: pmf_{i}+= p_{i}Δt