Issue when counting reaction in time instead of reaction steps

[lhuangb]

Hello!

I am trying to run a simulation where I am adding a high or low track at random, and I expect that for the same amount of time (such as t = 40), as my rate of high tracks increases, the number of high tracks will increase but the number of low tracks will stay the same. To implement this, I used an if t > 40 break statement, and when I run the function, it gives me the results as expected. However, when I apply the derivative estimate, when rate of high tracks increases, I am getting a negative value for rate of low tracks. This should be the result I get if I did not implement a restraint on time. Is my problem due to using an if statement? If so, is there another way I can ensure my reaction always stops at the specified time instead of number of reaction steps?
Any help is appreciated, thank you!

# Function
function Gillespie_trial3(rh1, rl1)
    
    s = [rand(1:1000) for i in 1:10000]
    rh = rh1 #rate of high tracks
    rl = rl1  #rate of low tracks
    Rtot = rh + rl
    t = 0 #s
    tracks = 0
    ltracks = 0
    htracks = 0
    times = Real[0]
    alltracks = Real[0]
    lowtracks = Real[0]
    hightracks = Real[0]
   
    for l in 1:500
        Random.seed!(s[l])
        dt = rand(Exponential(1/Rtot)) #time of next reaction
        probs = [(rh / (rh + rl)), (rl / (rh + rl))]
        Random.seed!(s[l])
        step_index = rand(Categorical(probs)) #high or low dose rate
        step = [1, 0][step_index]
        tracks += step #whole number of tracks
        t += dt
        
        if t > 40
            break
        end

        push!(times, t)
        push!(alltracks, tracks)
        htracks += step
        ltracks += 1-step
        push!(hightracks, htracks)
        push!(lowtracks, ltracks)
    end
    return htracks, ltracks, t
end

#Testing function to see results
x = [0.75, 1.5]
println(Gillespie_trial3(x[1], x[2]))
x = [1.5, 1.5]
println(Gillespie_trial3(x[1], x[2]))

#Derivative estimate
x = [0.75, 1.5]
g(p, q) = Gillespie_trial3(p, q)
function deriv2(x)
    a = derivative_estimate(l -> g(l, 1.5), x[1])
    b = derivative_estimate(m -> g(0.75, m), x[2])
    return a, b
end
deriv2(x)

[ChrisRackauckas]

If you add some logging on the t what do you see during the differentiation run?

[lhuangb]

Hello,
Thank you for your response! I tried adding logging to dt and stopped the reaction when t gets to a certain point (I chose -100 since my dts are small, so ln(dt) would be negative).

# Function
function Gillespie_trial3(rh1, rl1)
    
    s = [rand(1:1000) for i in 1:10000]
    rh = rh1 #rate of high tracks
    rl = rl1  #rate of low tracks
    Rtot = rh + rl
    t = 0 #s
    tracks = 0
    ltracks = 0
    htracks = 0
    times = Real[0]
    alltracks = Real[0]
    lowtracks = Real[0]
    hightracks = Real[0]
   
    for l in 1:500
        Random.seed!(s[l])
        dt = rand(Exponential(1/Rtot)) #time of next reaction
        probs = [(rh / (rh + rl)), (rl / (rh + rl))]
        Random.seed!(s[l])
        step_index = rand(Categorical(probs)) #high or low dose rate
        step = [1, 0][step_index]
        tracks += step #whole number of tracks
        t += log(dt)
        
        if t < -100
            break
        end

        push!(times, t)
        push!(alltracks, tracks)
        htracks += step
        ltracks += 1-step
        push!(hightracks, htracks)
        push!(lowtracks, ltracks)
    end
    return htracks, ltracks, t
end

#Testing function to see results
x = [0.75, 1.5]
println(Gillespie_trial3(x[1], x[2]))
x = [1.5, 1.5]
println(Gillespie_trial3(x[1], x[2]))

#Derivative estimate
x = [0.75, 1.5]
g(p, q) = Gillespie_trial3(p, q)
function deriv2(x)
    a = derivative_estimate(l -> g(l, 1.5), x[1])
    b = derivative_estimate(m -> g(0.75, m), x[2])
    return a, b
end
deriv2(x)

From this function, here is what I get:
For rate high = 0.75, rate low = 1.5:
(19, 51, -101.03453975640085) = (# high tracks, # low tracks, t)
For rate high = 1.5, rate low = 1.5
(28, 22, -102.14326202203881) = (# high tracks, # low tracks, t)

This is my differentiation when changing rate high while keeping rate low constant:
(16.000000000000004, -16.000000000000004, -27.555555555555518) = (differentiation for # high, # low, time)

I am getting the same differentiation pattern as when I did not use ln(dt), but it is not my desired differentiation.
For the function output, though the "t" values are similar, I wanted my # low tracks to stay around constant (shown below)

When I stop my reaction at t = 20 without log, here is what my function gives me:
For rate high = 0.75, rate low = 1.5:
(19, 33, 20.025476261383837) = (# high tracks, # low tracks, t)
For rate high = 1.5, rate low = 1.5:
(47, 35, 20.018093733837052) = (# high tracks, # low tracks, t)
Where the number of low tracks stays relatively constant while the number of high tracks has greatly increased.

When I do the differentiation run and change rate high while keeping rate low constant, here is what I get:
(13.777777777777782, -13.777777777777782, -9.08741795997318) = (differentiation for # high, # low, time)

So, it is a similar "pattern" to the differentiation with the log, but my function outputs are different.
As for my desired differentiation, I expect rate high to have a positive derivative, but I expect rate low to be closer to zero since the number of low tracks do not significantly change in my function output without log.
Therefore, I believe I still need a different method of stopping the differentiation at a fixed time so that my number of low tracks will stay relatively constant.
Please let me know if you have any suggestions, or if there is another way of adding logging to t that I am not thinking of! Thank you for your time.

[gaurav-arya]

Thanks for the issue! You are correct that the t > 40 check can result in incorrect derivative estimates, since t is a continuous random variable that depends on your parameters. See the fifth limitation on https://gaurav-arya.github.io/StochasticAD.jl/stable/limitations.html

One way to solve the issue is to use a bound on the total reaction rate that is independent of the parameters. We can then introduce a "ghost" reaction that has no effect (corresonding to step_index = 2 in the below code) to ensure that the rate used in the Gillespie algorithm is always equal to this bound. That way, although t is still random, it no longer depends on our parameters, so we get correct estimates. With the below implementation of this fix, the derivative estimates that we expect to be 0 are indeed exactly 0 all the time (although in general, we should only expect their expectation to be 0).

using StochasticAD
import Random
using Distributions

# Function
function Gillespie_trial3(rh1, rl1; Rtot_bound = rh1 + rl1)
    
    rh = rh1 #rate of high tracks
    rl = rl1  #rate of low tracks
    Rtot = rh + rl
    t = 0 #s
    tracks = 0
    ltracks = 0
    htracks = 0
   
    for l in 1:500
        dt = rand(Exponential(1/Rtot_bound)) #time of next reaction
        probs = [rh / Rtot_bound, 1 - Rtot / Rtot_bound, rl / Rtot_bound] 
        step_index = rand(Categorical(probs)) #high or low dose rate
        tracks += [1, 0, 1][step_index] #whole number of tracks
        t += dt
        
        if t > 40
            break
        end

        htracks += [1, 0, 0][step_index]
        ltracks += [0, 0, 1][step_index]
    end
    return htracks, ltracks, t
end

#Testing function to see results
x = [0.75, 1.5]
println(Gillespie_trial3(x[1], x[2]))
x = [1.5, 1.5]
println(Gillespie_trial3(x[1], x[2]))

#Derivative estimate
x = [0.75, 1.5]
function deriv2(x)
    a = derivative_estimate(l -> Gillespie_trial3(l, 1.5; Rtot_bound = sum(x) + 0.5), x[1])
    b = derivative_estimate(m -> Gillespie_trial3(0.75, m; Rtot_bound = sum(x) + 0.5), x[2])
    return a, b
end
mean(deriv2(x)[1][1] for i in 1:1000) # approximately 40
mean(deriv2(x)[1][2] for i in 1:1000) # 0
mean(deriv2(x)[2][1] for i in 1:1000) # 0
mean(deriv2(x)[2][2] for i in 1:1000) # approximately 40

[lhuangb]

This was a great suggestion and has been solving our problem so far, thank you very much!