SciML/JumpProcesses.jl 教程：离散随机过程与Gillespie方法详解

引言

本文将深入讲解如何使用SciML生态系统中的JumpProcesses.jl包来定义和模拟连续时间跳跃过程（也称为点过程）。在生物学领域，这类模型常被称为随机化学动力学模型或Gillespie模型。我们将通过一个经典的SIR传播模型示例，全面介绍JumpProcesses.jl的核心功能和使用方法。

模型背景：SIR传播动力学

SIR模型是流行病学中描述传播动态的经典模型，它将人群分为三类：

• 易感者(Susceptible, S)：可能被影响的健康人群
• 感染者(Infected, I)：已影响并可传播的群体
• 康复者(Recovered, R)：已恢复并获得免疫力的群体

模型包含两个基本反应过程：

1. 传播过程：易感者与感染者接触后转变为感染者
2. 康复过程：感染者康复并获得免疫力

数学建模

在数学上，SIR模型可以用跳跃过程来描述。设S(t)、I(t)、R(t)分别表示t时刻三类人群的数量，模型可表示为：

\begin{aligned}
S(t) &= S(0) - Y_1\left(\int_0^t \beta S(s^-)I(s^-)ds\right) \\
I(t) &= I(0) + Y_1\left(\int_0^t \beta S(s^-)I(s^-)ds\right) - Y_2\left(\int_0^t \nu I(s^-)ds\right) \\
R(t) &= R(0) + Y_2\left(\int_0^t \nu I(s^-)ds\right)
\end{aligned}

其中Y₁和Y₂是独立的单位泊松过程，β是传播率参数，ν是康复率参数。

实现方法一：使用Catalyst构建反应网络

对于化学/生物系统，推荐使用Catalyst.jl来方便地构建模型：

using Catalyst
sir_model = @reaction_network begin
    β, S + I --> 2I  # 传播过程
    ν, I --> R       # 康复过程
end

定义参数和初始条件：

p = (:β => 0.1/1000, :ν => 0.01)
u₀ = [:S => 990, :I => 10, :R => 0]
tspan = (0.0, 250.0)

构建离散问题并转换为跳跃问题：

prob = DiscreteProblem(sir_model, u₀, tspan, p)
jump_prob = JumpProblem(sir_model, prob, Direct())

求解并可视化：

sol = solve(jump_prob, SSAStepper())
plot(sol)

实现方法二：直接使用JumpProcesses底层接口

JumpProcesses提供了三种跳跃类型，按性能从高到低排序：

1. MassActionJump：性能最佳，适用于质量作用定律的反应
2. ConstantRateJump：中等性能，适用于速率恒定的跳跃
3. VariableRateJump：最灵活但性能最低，适用于速率随时间变化的跳跃

使用ConstantRateJump

定义传播和康复两个跳跃过程：

β = 0.1/1000.0
ν = 0.01
p = (β, ν)

# 传播过程
rate1(u, p, t) = p[1] * u[1] * u[2]  # β*S*I
function affect1!(integrator)
    integrator.u[1] -= 1  # S减少1
    integrator.u[2] += 1  # I增加1
end
jump1 = ConstantRateJump(rate1, affect1!)

# 康复过程
rate2(u, p, t) = p[2] * u[2]  # ν*I
function affect2!(integrator)
    integrator.u[2] -= 1  # I减少1
    integrator.u[3] += 1  # R增加1
end
jump2 = ConstantRateJump(rate2, affect2!)

构建并求解问题：

u₀ = [999, 10, 0]
tspan = (0.0, 250.0)
prob = DiscreteProblem(u₀, tspan, p)
jump_prob = JumpProblem(prob, Direct(), jump1, jump2)
sol = solve(jump_prob, SSAStepper())
plot(sol, label=["S(t)" "I(t)" "R(t)"])

性能优化建议

1. 对于符合质量作用定律的反应，优先使用MassActionJump
2. 当跳跃速率在跳跃间恒定不变时，使用ConstantRateJump
3. 只有当速率显式依赖时间时才考虑使用VariableRateJump
4. 对于VariableRateJump，如果可以提供速率上界，使用有界版本提高性能

扩展应用

JumpProcesses.jl还支持：

• 混合跳跃过程和连续动态（ODE/SDE）的混合系统
• 近似模拟的τ-leaping方法（通过RegularJump）
• 复杂的时间依赖速率函数
• 各种保存控制以减少内存使用

结论

本文详细介绍了使用JumpProcesses.jl构建和模拟离散随机过程的两种主要方法：通过Catalyst.jl的高级接口和直接使用JumpProcesses的底层接口。通过SIR传播模型的示例，我们展示了如何定义跳跃过程、设置模拟参数以及可视化结果。理解不同跳跃类型的性能特点和适用场景，可以帮助用户在模型复杂性和计算效率之间做出合理的选择。

JumpProcesses.jl作为SciML生态系统的重要组成部分，为复杂随机系统的建模和仿真提供了强大而灵活的工具，特别适合生物化学、流行病学、系统生物学等领域的研究应用。