Python实现季节性驱动SEIR模型分岔图及教材复现疑问

如何在Python中实现SEIR（易感、暴露、感染、恢复）这类季节性驱动流行病学模型的分岔图？已掌握该模型的实现及时间序列采样展示方法，但难以复现Keeling 2007年著作中的分岔图。

正在复现《Modeling Infectious Diseases in Humans and Animals》中的图表，验证模型正确性并可视化参数对动力系统的影响。参考logistic映射分岔图实现，总结出分岔图的x分量为参数取值、y分量为对应时间点的模型解，据此实现了plot_bifurcation函数，但面板输出与教材图表不符，分岔图结果异常，代码执行无警告或错误。

教材分岔图：

我的面板输出：

我的分岔图输出：

1. 为何面板输出与教材图表不符？
2. 分岔图函数输出结果异常，是否对分岔图的理解存在误区？

1. 参数单位错误（核心问题） #

模型时间步为天，但设置的beta_zero=1250是年度传输率，未转换为天单位，导致实际传输率被放大365倍，完全偏离教材参数。需将beta_zero修正为1250 / days_per_year，确保所有速率单位统一为天。

2. 分岔图绘制逻辑误区 #

分岔图不应绘制所有时间点的解，而是要：

• 先让系统运行足够长时间，收敛到稳态/周期态
• 提取收敛后周期内的极值点（峰值/谷值），而非所有时间点的数值
• 对每个参数，仅绘制这些极值点，才能呈现分岔结构

3. 系统收敛时间不足 #

需先运行足够长的“预热期”（比如100年），让系统从初始条件收敛到稳定周期或混沌态，再采样后续时间序列的极值，避免初始暂态影响结果。

from typing import Callable, Dict, List, Optional, Any
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def seasonal_seir(y: List, t: List, params: Dict[str, Any]):
    """Seasonally forced SEIR model.
    
    Function parameters match those required by `scipy.integrate.odeint`

    Args:
        y: Initial conditions.
        t: Timesteps over which numerical solution will be computed.
        params: Dict with key-value pairs:
            beta_zero -- Average daily transmission rate.
            beta_one  -- Amplitude of seasonal forcing.
            omega     -- Period of forcing (2pi/365 for annual).
            mu        -- Daily natural mortality rate.
            sigma     -- Daily latent period inverse (1/latent_days).
            gamma     -- Daily recovery rate inverse (1/recovery_days).

    Returns:
        Tuple of derivatives for susceptible, exposed, infected variables.
    """
    beta_zero = params['beta_zero']
    beta_one = params['beta_one']
    omega = params['omega']
    mu = params['mu']
    sigma = params['sigma']
    gamma = params['gamma']

    s, e, i = y 
    beta = beta_zero*(1 + beta_one*np.cos(omega*t))
    sdot = mu - (beta * i + mu)*s
    edot = beta*s*i - (mu + sigma)*e
    idot = sigma*e - (mu + gamma)*i
    return sdot, edot, idot

def plot_panels(
    model: Callable,
    model_params: Dict, 
    panel_param_space: List, 
    panel_param_name: str,
    initial_conditions: List,
    warmup_timesteps: List,
    sample_timesteps: List,
    odeint_kwargs: Optional[Dict] = dict(),
    state_var_ix: Optional[int] = None,
    log_scale: bool = False):
    """Plot time series panels after system convergence."""

    fig, axs = plt.subplots(ncols=len(panel_param_space), figsize=(15,8))
    sample_years = sample_timesteps / days_per_year

    for ix, panel_param in enumerate(panel_param_space):
        model_params[panel_param_name] = panel_param

        # 先运行预热期，让系统收敛
        warmup_sol = odeint(
            model,
            initial_conditions,
            warmup_timesteps,
            args=(model_params,),
            **odeint_kwargs)
        # 用预热期的最终状态作为采样期的初始条件
        final_warmup_state = warmup_sol[-1]
        sample_sol = odeint(
            model,
            final_warmup_state,
            sample_timesteps,
            args=(model_params,),
            **odeint_kwargs)

        if state_var_ix is not None:
            sample_sol = sample_sol[:, state_var_ix]
        else:
            raise ValueError("Must specify state_var_ix for multi-variable model")

        if log_scale:
            sample_sol = np.log(sample_sol)

        axs[ix].plot(sample_years, sample_sol)
        axs[ix].set_title(f"{panel_param_name} = {panel_param}")

    fig.suptitle('Reproduced Panels from Keeling 2007')
    fig.supxlabel('Time (years)')
    fig.supylabel('Fraction Infected')
    plt.tight_layout()
    return fig, axs

def plot_bifurcation(
    model: Callable,
    model_params: Dict,
    bifurcation_parameter_space: np.ndarray,
    bifurcation_param_name: str,
    initial_conditions: List,
    warmup_timesteps: List,
    sample_timesteps: List,
    odeint_kwargs: Optional[Dict] = dict(),
    state_var_ix: Optional[int] = None,
    log_scale: bool = False):
    """Plot bifurcation diagram using extreme values after convergence."""

    fig, ax = plt.subplots(figsize=(15,8))

    for param in bifurcation_parameter_space:
        model_params[bifurcation_param_name] = param

        # 预热期收敛
        warmup_sol = odeint(
            model,
            initial_conditions,
            warmup_timesteps,
            args=(model_params,),
            **odeint_kwargs)
        final_warmup_state = warmup_sol[-1]
        sample_sol = odeint(
            model,
            final_warmup_state,
            sample_timesteps,
            args=(model_params,),
            **odeint_kwargs)

        if state_var_ix is not None:
            sample_sol = sample_sol[:, state_var_ix]
        else:
            raise ValueError("Must specify state_var_ix for multi-variable model")

        if log_scale:
            sample_sol = np.log(sample_sol)

        # 提取每年的峰值（假设采样期是整数年）
        annual_peaks = []
        n_sample_years = len(sample_timesteps) // days_per_year
        for year in range(n_sample_years):
            start_idx = year * days_per_year
            end_idx = (year+1)*days_per_year
            annual_segment = sample_sol[start_idx:end_idx]
            annual_peaks.append(np.max(annual_segment))
            # 可选：添加谷值
            # annual_peaks.append(np.min(annual_segment))

        # 绘制当前参数对应的所有极值点
        ax.scatter([param]*len(annual_peaks), annual_peaks, s=1, c='k', alpha=0.5)

    ax.set_title('Reproduced Bifurcation Diagram from Keeling 2007')
    ax.set_xlabel(r'$\beta_1$')
    ax.set_ylabel('Fraction Infected (Peak Values)')
    plt.tight_layout()
    return fig, ax

# ----------------------参数设置----------------------
days_per_year = 365

# 修正参数：所有速率统一为天单位
mu = 0.02 / days_per_year          # 日死亡率（年度0.02）
beta_zero = 1250 / days_per_year   # 日传输率（年度1250）
sigma = 1/8                        # 日潜伏率（潜伏期8天）
gamma = 1/5                        # 日恢复率（恢复时间5天）
omega = 2*np.pi / days_per_year    # 年度周期角频率

model_params = dict(
        beta_zero=beta_zero,
        omega=omega,
        mu=mu,
        sigma=sigma,
        gamma=gamma)

# 初始条件
s0 = 6e-2
e0 = i0 = 1e-3
initial_conditions = [s0, e0, i0]

# 时间步设置：100年预热期，100年采样期
warmup_nyears = 100
warmup_timesteps = np.arange(0, warmup_nyears*days_per_year, 1)

sample_nyears = 100
sample_timesteps = np.arange(0, sample_nyears*days_per_year, 1)

# 面板参数
beta_ones = [0.025, 0.05, 0.25]
infection_ix = 2  # 感染状态变量索引

# 绘制面板图
panel_fig, panel_ax = plot_panels(
    model=seasonal_seir,
    model_params=model_params,
    panel_param_space=beta_ones,
    panel_param_name='beta_one',
    initial_conditions=initial_conditions,
    warmup_timesteps=warmup_timesteps,
    sample_timesteps=sample_timesteps,
    odeint_kwargs=dict(hmax=5),
    state_var_ix=infection_ix,
    log_scale=False)

# 分岔图参数空间
bifurcation_beta_ones = np.linspace(0.025, 0.3, 200)

# 绘制分岔图
bi_fig, bi_ax = plot_bifurcation(
    model=seasonal_seir,
    model_params=model_params,
    bifurcation_parameter_space=bifurcation_beta_ones,
    bifurcation_param_name='beta_one',
    initial_conditions=initial_conditions,
    warmup_timesteps=warmup_timesteps,
    sample_timesteps=sample_timesteps,
    odeint_kwargs={'hmax':5},
    state_var_ix=infection_ix,
    log_scale=False)

plt.show()

• 参数修正后，面板图将呈现教材中的周期特征：β₁=0.025为年度单周期，β₁=0.05为2年周期，β₁=0.25为混沌态
• 分岔图通过提取每年峰值绘制，将清晰展示从周期到混沌的分岔过程
• 预热期确保系统收敛到稳定态，避免初始条件的暂态干扰

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