mutable struct JumpProblem{iip, P, A, C, J<:Union{Nothing, JumpProcesses.AbstractJumpAggregator}, J2, J3, J4, R, K} <: SciMLBase.AbstractJumpProblem{P, J<:Union{Nothing, JumpProcesses.AbstractJumpAggregator}}

Defines a collection of jump processes to associate with another problem type.

• Documentation Page
• Tutorial Page
• FAQ Page

Constructors

JumpProblems can be constructed by first building another problem type to which the jumps will be associated. For example, to simulate a collection of jump processes for which the transition rates are constant between jumps (called ConstantRateJumps or MassActionJumps), we must first construct a DiscreteProblem

prob = DiscreteProblem(u0, p, tspan)

where u0 is the initial condition, p the parameters and tspan the time span. If we wanted to have the jumps coupled with a system of ODEs, or have transition rates with explicit time dependence, we would use an ODEProblem instead that defines the ODE portion of the dynamics.

Given prob we define the jumps via

• JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::JumpSet ; kwargs...)
• JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps...; kwargs...)

Here aggregator specifies the underlying algorithm for calculating next jump times and types, for example Direct. The collection of different AbstractJump types can then be passed within a single JumpSet or as subsequent sequential arguments.

Fields

• prob: The type of problem to couple the jumps to. For a pure jump process use DiscreteProblem, to couple to ODEs, ODEProblem, etc.

• aggregator: The aggregator algorithm that determines the next jump times and types for ConstantRateJumps and MassActionJumps. Examples include Direct.

• discrete_jump_aggregation: The underlying state data associated with the chosen aggregator.

• jump_callback: CallBackSet with the underlying ConstantRate and VariableRate jumps.

• variable_jumps: The VariableRateJumps.

• regular_jump: The RegularJumps.

• massaction_jump: The MassActionJumps.

• rng: The random number generator to use.

• kwargs: kwargs to pass on to solve call.

Keyword Arguments

• rng, the random number generator to use. Defaults to Julia's built-in generator.
• save_positions=(true,true) when including variable rates and (false,true) for constant rates, specifies whether to save the system's state (before, after) the jump occurs.
• spatial_system, for spatial problems the underlying spatial structure.
• hopping_constants, for spatial problems the spatial transition rate coefficients.
• use_vrj_bounds = true, set to false to disable handling bounded VariableRateJumps with a supporting aggregator (such as Coevolve). They will then be handled via the continuous integration interface, and treated like general VariableRateJumps.
• vr_aggregator, indicates the aggregator to use for sampling variable rate jumps. Current default is VR_FRM.

Please see the tutorial page in the DifferentialEquations.jl docs for usage examples and commonly asked questions.

struct SSAStepper <: SciMLBase.AbstractDEAlgorithm

Highly efficient integrator for pure jump problems that involve only ConstantRateJumps, MassActionJumps, and/or VariableRateJumps with rate bounds.

Notes

• Only works with JumpProblems defined from DiscreteProblems.
• Only works with collections of ConstantRateJumps, MassActionJumps, and VariableRateJumps with rate bounds.
• Only supports DiscreteCallbacks for events, which are checked after every step taken by SSAStepper.
• Only supports a limited subset of the output controls from the common solver interface, specifically save_start, save_end, and saveat.
• As when using jumps with ODEs and SDEs, saving controls for whether to save each time a jump occurs are via the save_positions keyword argument to JumpProblem. Note that when choosing SSAStepper as the timestepper, save_positions = (true,true), (true,false), or (false,true) are all equivalent. SSAStepper will save only the post-jump state in the solution object in each of these cases. This is because solution objects generated via SSAStepper use piecewise-constant interpolation, and can therefore exactly reconstruct the sampled jump process path with knowing just the post-jump state. That is, sol(t) for any 0 <= t <= tstop will give the exact value of the sampled solution path at t provided at least one component of save_positions is true.

Examples

SIR model:

using JumpProcesses
β = 0.1 / 1000.0; ν = .01;
p = (β,ν)
rate1(u,p,t) = p[1]*u[1]*u[2]  # β*S*I
function affect1!(integrator)
  integrator.u[1] -= 1         # S -> S - 1
  integrator.u[2] += 1         # I -> I + 1
end
jump = ConstantRateJump(rate1,affect1!)

rate2(u,p,t) = p[2]*u[2]      # ν*I
function affect2!(integrator)
  integrator.u[2] -= 1        # I -> I - 1
  integrator.u[3] += 1        # R -> R + 1
end
jump2 = ConstantRateJump(rate2,affect2!)
u₀    = [999,1,0]
tspan = (0.0,250.0)
prob = DiscreteProblem(u₀, tspan, p)
jump_prob = JumpProblem(prob, Direct(), jump, jump2)
sol = solve(jump_prob, SSAStepper())

see the tutorial for details.

reset_aggregated_jumps!(integrator, uprev = nothing; update_jump_params=true)

Reset the state of jump processes and associated solvers following a change in parameters or such.

Notes

• update_jump_params=true will recalculate the rates stored within any MassActionJump that was built from the parameter vector. If the parameter vector is unchanged, this can safely be set to false to improve performance.

struct ConstantRateJump{F1, F2} <: JumpProcesses.AbstractJump

Defines a jump process with a rate (i.e. hazard, intensity, or propensity) that does not explicitly depend on time. More precisely, one where the rate function is constant between the occurrence of jumps. For detailed examples and usage information, see the

• Tutorial

Fields

• rate: Function rate(u,p,t) that returns the jump's current rate.

• affect!: Function affect(integrator) that updates the state for one occurrence of the jump.

Examples

Suppose u[1] gives the number of particles and p[1] the probability per time each particle can decay away. A corresponding ConstantRateJump for this jump process is

rate(u,p,t) = p[1]*u[1]
affect!(integrator) = integrator.u[1] -= 1
crj = ConstantRateJump(rate, affect!)

Notice, here that rate changes in time, but is constant between the occurrence of jumps (when u[1] will decrease).

struct MassActionJump{T, S, U, V} <: JumpProcesses.AbstractMassActionJump

Optimized representation for ConstantRateJumps that can be represented in mass action form, offering improved performance within jump algorithms compared to ConstantRateJump. For detailed examples and usage information, see the

• Main Docs
• Tutorial

Constructors

• MassActionJump(reactant_stoich, net_stoich; scale_rates = true, param_idxs = nothing)

Here reactant_stoich denotes the reactant stoichiometry for each reaction and net_stoich the net stoichiometry for each reaction.

Fields

• scaled_rates: The (scaled) reaction rate constants.

• reactant_stoch: The reactant stoichiometry vectors.

• net_stoch: The net stoichiometry vectors.

• param_mapper: Parameter mapping functor to identify reaction rate constants with parameters in p vectors.

Keyword Arguments

• scale_rates = true, whether to rescale the reaction rate constants according to the stoichiometry.
• nocopy = false, whether the MassActionJump can alias the scaled_rates and reactant_stoch from the input. Note, if scale_rates=true this will potentially modify both of these.
• param_idxs = nothing, indexes in the parameter vector, JumpProblem.prob.p, that correspond to each reaction's rate.

See the tutorial and main docs for details.

Examples

An SIR model with S + I --> 2I at rate β as the first reaction and I --> R at rate ν as the second reaction can be encoded by

p        = (β=1e-4, ν=.01)
u0       = [999, 1, 0]       # (S,I,R)
tspan    = (0.0, 250.0)
rateidxs = [1, 2]           # i.e. [β,ν]
reactant_stoich =
[
  [1 => 1, 2 => 1],         # 1*S and 1*I
  [2 => 1]                  # 1*I
]
net_stoich =
[
  [1 => -1, 2 => 1],        # -1*S and 1*I
  [2 => -1, 3 => 1]         # -1*I and 1*R
]
maj = MassActionJump(reactant_stoich, net_stoich; param_idxs=rateidxs)
prob = DiscreteProblem(u0, tspan, p)
jprob = JumpProblem(prob, Direct(), maj)

Notes

• By default, reaction rates are rescaled when constructing the MassActionJump as explained in the main docs. Disable this with the kwarg scale_rates=false.
• Also see the main docs for how to specify reactions with no products or no reactants.

struct VariableRateJump{R, F, R2, R3, R4, I, T, T2} <: JumpProcesses.AbstractJump

Defines a jump process with a rate (i.e. hazard, intensity, or propensity) that may explicitly depend on time. More precisely, one where the rate function is allowed to change between the occurrence of jumps. For detailed examples and usage information, see the

• Tutorial

Note that two types of VariableRateJumps are currently supported, with different performance characteritistics.

• A general VariableRateJump or VariableRateJump will refer to one in which only rate and affect functions are specified.

  • These are the most general in what they can represent, but require the use of an ODEProblem or SDEProblem whose underlying timestepper handles their evolution in time (via the callback interface).
  • This is the least performant jump type in simulations.

• Bounded VariableRateJumps require passing the keyword arguments urate and rateinterval, corresponding to functions urate(u, p, t) and rateinterval(u, p, t), see below. These must calculate a time window over which the rate function is bounded by a constant. Note that it is ok if the rate bound would be violated within the time interval due to a change in u arising from another ConstantRateJump, MassActionJump or boundedVariableRateJump being executed, as the chosen aggregator will then handle recalculating the rate bound and interval. However, if the bound could be violated within the time interval due to a change in u arising from continuous dynamics such as a coupled ODE, SDE, or a general VariableRateJump, bounds should not be given. This ensures the jump is classified as a general VariableRateJump and properly handled. One can also optionally provide a lower bound function, lrate(u, p, t), via the lrate keyword argument. This can lead to increased performance. The validity of the lower bound should hold under the same conditions and rate interval as urate.

  • Bounded VariableRateJumps can currently be used in the Coevolve aggregator, and can therefore be efficiently simulated in pure-jump DiscreteProblems using the SSAStepper time-stepper.
  • These can be substantially more performant than general VariableRateJumps without the rate bound functions.

Reemphasizing, the additional user provided functions leveraged by bounded VariableRateJumps, urate(u, p, t), rateinterval(u, p, t), and the optional lrate(u, p, t) require that

• For s in [t, t + rateinterval(u, p, t)], we have that lrate(u, p, t) <= rate(u, p, s) <= urate(u, p, t).
• It is ok if these bounds would be violated during the time window due to another ConstantRateJump, MassActionJump or bounded VariableRateJump occurring. However, they must remain valid if u changes for any other reason (for example, due to continuous dynamics like ODEs, SDEs, or general VariableRateJumps).

Fields

• rate: Function rate(u,p,t) that returns the jump's current rate given state u, parameters p and time t.

• affect!: Function affect!(integrator) that updates the state for one occurrence of the jump given integrator.

• lrate: Optional function lrate(u, p, t) that computes a lower bound on the rate in the interval t to t + rateinterval(u, p, t) at time t given state u and parameters p. This bound must rigorously hold during the time interval as long as another ConstantRateJump, MassActionJump, or boundedVariableRateJump has not been sampled. When using aggregators that support bounded VariableRateJumps, currently only Coevolve, providing a lower-bound can lead to improved performance.

• urate: Optional function urate(u, p, t) for general VariableRateJumps, but is required to define a bounded VariableRateJump, which can be used with supporting aggregators, currently only Coevolve, and offers improved computational performance. Computes an upper bound for the rate in the interval t to t + rateinterval(u, p, t) at time t given state u and parameters p. This bound must rigorously hold during the time interval as long as another ConstantRateJump, MassActionJump, or boundedVariableRateJump has not been sampled.

• rateinterval: Optional function rateinterval(u, p, t) for general VariableRateJumps, but is required to define a bounded VariableRateJump, which can be used with supporting aggregators, currently only Coevolve, and offers improved computational performance. Computes the time interval from time t over which the urate and lrate bounds will hold, t to t + rateinterval(u, p, t), given state u and parameters p. This bound must rigorously hold during the time interval as long as another ConstantRateJump, MassActionJump, or boundedVariableRateJump has not been sampled.

• idxs

• rootfind

• interp_points

• save_positions

• abstol

• reltol

Examples

Suppose u[1] gives the number of particles and t*p[1] the probability per time each particle can decay away. A corresponding VariableRateJump for this jump process is

rate(u,p,t) = t*p[1]*u[1]
affect!(integrator) = integrator.u[1] -= 1
vrj = VariableRateJump(rate, affect!)

To define a bounded VariableRateJump that can be used with supporting aggregators such as Coevolve, we must define bounds and a rate interval:

rateinterval(u,p,t) = (1 / p[1]) * 2
rate(u,p,t) = t * p[1] * u[1]
lrate(u, p, t) = rate(u, p, t)
urate(u,p,t) = rate(u, p, t + rateinterval(u,p,t))
affect!(integrator) = integrator.u[1] -= 1
vrj = VariableRateJump(rate, affect!; lrate = lrate, urate = urate,
                                      rateinterval = rateinterval)

Notes

• When using an aggregator that supports bounded VariableRateJumps, DiscreteProblem can be used. Otherwise, ODEProblem or SDEProblem must be used.
• When not using aggregators that support bounded VariableRateJumps, or when there are general VariableRateJumps, integrators store an effective state type that wraps the main state vector. See ExtendedJumpArray for details on using this object. In this case all ConstantRateJump, VariableRateJump and callback affect! functions receive an integrator with integrator.u an ExtendedJumpArray.
• Salis H., Kaznessis Y., Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions, Journal of Chemical Physics, 122 (5), DOI:10.1063/1.1835951 is used for calculating jump times with VariableRateJumps within ODE/SDE integrators.

struct RegularJump{iip, R, C, MD}

Representation for encoding rates and multiple simultaneous jumps for use in τ-leaping type methods.

Constructors

• RegularJump(rate, c, numjumps; mark_dist = nothing)

Fields

• rate: Function rate!(rate_vals, u, p, t) that returns the current rates, i.e. intensities or propensities, for all possible jumps in rate_vals.

• c: Function c(du, u, p, t, counts, mark) that executes the ith jump counts[i] times, saving the output in du[i].

• numjumps: Number of jumps in the system.

• mark_dist: A distribution for marks. Not currently used or supported.

Examples

function rate!(out, u, p, t)
    out[1] = (0.1 / 1000.0) * u[1] * u[2]
    out[2] = 0.01u[2]
    nothing
end

const dc = zeros(3,2)
function c(du, u, p, t, counts, mark) 
    mul!(du, dc, counts)
    nothing
end

rj = RegularJump(rate!, c, 2)

## Notes
- `mark_dist` is not currently used or supported in τ-leaping methods.

struct JumpSet{T1, T2, T3, T4} <: JumpProcesses.AbstractJump

Defines a collection of jumps that should collectively be included in a simulation.

Fields

• variable_jumps: Collection of VariableRateJumps

• constant_jumps: Collection of ConstantRateJumps

• regular_jump: Collection of RegularJumps

• massaction_jump: Collection of MassActionJumps

Examples

Here we construct two jumps, store them in a JumpSet, and then simulate the resulting process.

using JumpProcesses, OrdinaryDiffEq

rate1(u,p,t) = p[1]
affect1!(integrator) = (integrator.u[1] += 1)
crj = ConstantRateJump(rate1, affect1!)

rate2(u,p,t) = (t/(1+t))*p[2]*u[1]
affect2!(integrator) = (integrator.u[1] -= 1)
vrj = VariableRateJump(rate2, affect2!)

jset = JumpSet(crj, vrj)

f!(du,u,p,t) = (du .= 0)
u0 = [0.0]
p = (20.0, 2.0)
tspan = (0.0, 200.0)
oprob = ODEProblem(f!, u0, tspan, p)
jprob = JumpProblem(oprob, Direct(), jset)
sol = solve(jprob, Tsit5())

An improvement of the COEVOLVE algorithm for simulating any compound jump process that evolves through time. This method handles variable intensity rates with user-defined bounds and inter-dependent processes. It reduces to NRM when rates are constant. As opposed to COEVOLVE, this method syncs the thinning procedure with the stepper which allows it to handle dependencies on continuous dynamics.

G. A. Zagatti, S. A. Isaacson, C. Rackauckas, V. Ilin, S.-K. Ng and S. Bressan, Extending JumpProcess.jl for fast point process simulation with time-varying intensities, arXiv. doi:10.48550/arXiv.2306.06992.

M. Farajtabar, Y. Wang, M. Gomez-Rodriguez, S. Li, H. Zha, and L. Song, COEVOLVE: a joint point process model for information diffusion and network evolution, Journal of Machine Learning Research 18(1), 1305–1353 (2017). doi: 10.5555/3122009.3122050.

Gillespie's Direct method. ConstantRateJump rates and affects are stored in tuples. Fastest for a small (total) number of ConstantRateJumps or MassActionJumps (~10). For larger numbers of possible jumps, use other methods.

Daniel T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, Journal of Computational Physics, 22 (4), 403–434 (1976). doi:10.1016/0021-9991(76)90041-3.

The Composition-Rejection Direct method. Performs best relative to other methods for systems with large numbers of jumps with special structure (for example a linear chain of reactions, or jumps corresponding to particles hopping on a grid or graph).

A. Slepoy, A.P. Thompson and S.J. Plimpton, A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks, Journal of Chemical Physics, 128 (20), 205101 (2008). doi:10.1063/1.2919546.

S. Mauch and M. Stalzer, Efficient formulations for exact stochastic simulation of chemical systems, ACM Transactions on Computational Biology and Bioinformatics, 8 (1), 27-35 (2010). doi:10.1109/TCBB.2009.47.

Gillespie's First Reaction Method. Should not be used for practical applications due to slow performance relative to all other methods.

Daniel T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, Journal of Computational Physics, 22 (4), 403–434 (1976). doi:10.1016/0021-9991(76)90041-3.

The Next Reaction Method. Can significantly outperform Direct for systems with large numbers of jumps and sparse dependency graphs, but is usually slower than one of DirectCR, RSSA, or RSSACR for such systems.

M. A. Gibson and J. Bruck, Efficient exact stochastic simulation of chemical systems with many species and many channels, Journal of Physical Chemistry A, 104 (9), 1876-1889 (2000). doi:10.1021/jp993732q.

A rejection-based direct method.

The Rejection SSA method. One of the best methods for systems with hundreds to many thousands of jumps (along with RSSACR) and sparse dependency graphs.

V. H. Thanh, C. Priami and R. Zunino, Efficient rejection-based simulation of biochemical reactions with stochastic noise and delays, Journal of Chemical Physics, 141 (13), 134116 (2014). doi:10.1063/1.4896985

V. H. Thanh, R. Zunino and C. Priami, On the rejection-based algorithm for simulation and analysis of large-scale reaction networks, Journal of Chemical Physics, 142 (24), 244106 (2015). doi:10.1063/1.4922923.

The Rejection SSA Composition-Rejection method. Often the best performer for systems with tens of thousands of jumps and sparse dependency graphs.

V. H. Thanh, R. Zunino, and C. Priami, Efficient constant-time complexity algorithm for stochastic simulation of large reaction networks, IEEE/ACM Transactions on Computational Biology and Bioinformatics, 14 (3), 657-667 (2017). doi:10.1109/TCBB.2016.2530066.

The Sorting Direct method. Often the fastest algorithm for smaller to moderate sized systems (tens of jumps), or systems where a few jumps occur much more frequently than others.

J. M. McCollum, G. D. Peterson, C. D. Cox, M. L. Simpson and N. F. Samatova, The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior, Computational Biology and Chemistry, 30 (1), 39049 (2006). doi:10.1016/j.compbiolchem.2005.10.007.

struct ExtendedJumpArray{T3<:Number, T1, T<:AbstractArray{T3<:Number, T1}, T2} <: AbstractArray{T3<:Number, 1}

Extended state definition used within integrators when there are VariableRateJumps in a system. For detailed examples and usage information, see the

• Tutorial

Fields

• u: The current state.

• jump_u: The current rate (i.e. hazard, intensity, or propensity) values for the VariableRateJumps.

Examples

using JumpProcesses, OrdinaryDiffEq
f(du,u,p,t) = du .= 0
rate(u,p,t) = (1+t)*u[1]*u[2]

# suppose we wish to decrease each of the two variables by one
# when a jump occurs
function affect!(integrator)
   # Method 1, direct indexing works like normal
   integrator.u[1] -= 1
   integrator.u[2] -= 1

   # Method 2, if we want to broadcast or use array operations we need
   # to access integrator.u.u which is the actual state object.
   # So equivalently to above we could have said:
   # integrator.u.u .-= 1
end

u0 = [10.0, 10.0]
vrj = VariableRateJump(rate, affect!)
oprob = ODEProblem(f, u0, (0.0,2.0))
jprob = JumpProblem(oprob, Direct(), vrj)
sol = solve(jprob,Tsit5())

Notes

• If ueja isa ExtendedJumpArray with ueja.u of size N and ueja.jump_u of size num_variableratejumps then

  # for 1 <= i <= N
  ueja[i] == ueja.u[i]
  
  # for N < i <= (N+num_variableratejumps)
  ueja[i] == ueja.jump_u[i]

• In a system with VariableRateJumps all callback, ConstantRateJump, and VariableRateJumpaffect! functions will receive integrators with integrator.u an ExtendedJumpArray.

• As such, affect! functions that wish to modify the state via vector operations should use ueja.u.u to obtain the aliased state object.

mutable struct SSAIntegrator{F, uType, tType, tdirType, P, S, CB, SA, OPT, TS} <: JumpProcesses.AbstractSSAIntegrator{SSAStepper, Nothing, uType, tType}

Solution objects for pure jump problems solved via SSAStepper.

Fields

• f: The underlying prob.f function. Not currently used.

• u: The current solution values.

• t: The current solution time.

• tprev: The previous time a jump occurred.

• tdir: The direction time is changing in (must be positive, indicating time is increasing)

• p: The current parameters.

• sol: The current solution object.

• i

• tstop: The next jump time.

• cb: The jump aggregator callback.

• saveat: Times to save the solution at.

• save_everystep: Whether to save every time a jump occurs.

• save_end: Whether to save at the final step.

• cur_saveat: Index of the next saveat time.

• opts: Tuple storing callbacks.

• tstops: User supplied times to step to, useful with callbacks.

• tstops_idx

• u_modified

• keep_stepping

• alias_tstops: If true, will write tstops into the user-passed array

• copied_tstops: If true indicates we have already allocated the tstops array