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?