元种群SIR模型中易感者(S)出现负值的问题求助

我正在研究基于元种群的SIR（易感-感染-康复）模型，模型通过基于距离的评分实现细胞间两两关联。运行时发现易感者状态(S)出现负值，无法明确原因。

推测这可能是因为未对变量变化设置约束导致，想确认是否是自身代码逻辑问题，同时寻求解决方案。

我制作了一个简化的复现代码示例，用随机值模拟场景：当感染速率相对较小时模型运行正常，但当速率因大量邻近细胞被感染（此处用随机数模拟突变）时，求解器估计失败。

library(deSolve)
ode(
  y = c(
    S1 = 10,
    I1 = 0.001),
  times = 0:5000,
  parms = list(# beta = 0.05
    beta = 0.01),
  func = function(times, state, parms) {
    with(parms, {
      # S <- state[c(1, 3)]
      # I <- state[c(2, 4)]
      S <- state[1]
      I <- state[2]

      beta <- beta + 2 * rbinom(1, size = 1, prob = 0.4)

      list(c(dS = -beta * S * I,
             dI = beta * S * I))
    })
  },

  jacfunc = function(t, y, p) {
    with(p, {
      S <- y[1]
      I <- y[2]
      (matrix(
        c(beta * I, beta * S, -beta * I, -beta * S),
        nrow = 2,
        ncol = 2,
        byrow = TRUE
      ))
    })
  },
  jactype = "fullusr",
  verbose = TRUE
) %>%
  # diagnostics() %>%
  print() %>%
  # zapsmall()
  identity()

已尝试的无效方法：

• 指定hmin参数无法解决该问题
• 提供自定义雅可比矩阵也无济于事

内容的提问来源于stack exchange，提问作者Mossa