SEIQR传染病模型Python实现行为异常问题排查与修正求助

各位好！我正在做微分方程课程的项目，需要把一个现有的SEIQR传染病传播模型用Python实现出来。我找到了一套带有真实参数的模型，但自己写的代码运行后，结果完全不对——最离谱的是，有些种群的曲线数值居然超过了总人口数，这显然违反了基本逻辑。我现在完全摸不着头脑，不知道是方程写错了、参数设错了，还是有其他问题，也不确定能不能在不破坏模型动力学特性的前提下，把数值限制在合理范围内。

下面是我的代码，麻烦大家帮忙看看哪里出问题了：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def seiqr_model(y, t, vinv, beta1, beta2, sigma1, sigma2, sigma3, r, d1, d2, alpha, eta1, eta3):
    S, E, I, R, Q = y
    dSdt = vinv - d1*S - alpha*vinv/d1 * E
    dEdt = ((alpha*vinv - d1 *eta1)/d1) * E
    dIdt = r*E - eta3*I
    dRdt = sigma3*E + sigma2*I - d1*R + sigma1*Q
    dQdt = beta1*E + beta2*I - eta3*Q
    return [dSdt, dEdt, dIdt, dRdt, dQdt]

# 初始条件
S0 = 34218169
E0 = 10000
I0 = 157
R0 = 0
Q0 = 0
y0 = [S0, E0, I0, R0, Q0,]

# 参数设置
vinv = 2300
beta1 = 0.15
beta2 = 0.00001
sigma1 = 0.00001
sigma2 = 0.00001
sigma3 = 0.001
r = 0.00325
d1 = 0.00003
d2 = 0.0000003423
alpha = 0.00000000264
eta1 = r+beta1+sigma3+d1
eta3 = sigma1+d1+d2

t = np.linspace(0, 75, 75)
sol = odeint(seiqr_model, y0, t, args=(vinv, beta1, beta2, sigma1, sigma2, sigma3, r, d1, d2, alpha, eta1, eta3))

plt.plot(t, sol[:,0], label='S')
plt.plot(t, sol[:,1], label='E')
plt.plot(t, sol[:,2], label='I')
plt.plot(t, sol[:,3], label='R')
plt.plot(t, sol[:,4], label='Q')
plt.xlabel('Time')
plt.ylabel('Proportion')
plt.legend()
plt.show()

我现在的几个疑问：

• 是不是微分方程的推导有错误？比如各种群的变化率公式哪里写错了，导致数值溢出？
• 如果要让所有种群数值不超过总人口，有没有办法在不改变模型核心动力学的前提下实现？
• 参数的设置是不是有问题？比如某些参数的量级或者取值逻辑不对，导致模型行为异常？

真心希望懂传染病模型或者微分方程数值求解的朋友能帮我分析一下，谢谢大家！

备注：内容来源于stack exchange，提问作者Ali Ismail