如何为基因切换开关模型绘制分岔图？求MATLAB实现指引

核心思路 #

对于2D切换开关系统，双稳态区域的边界对应鞍结分岔——此时系统的平衡点从「2个稳定点+1个鞍点」变为「1个稳定点」（或反之），雅可比矩阵会出现一个零特征值。我们需要遍历α₁和α₂的取值范围，对每个参数组合找到所有平衡点，通过雅可比矩阵判断平衡点稳定性，最终确定双稳态区域的边界。

MATLAB实现步骤提示 #

1. 定义参数遍历范围 #

先确定α₁和α₂的区间（示例取0到10，步长0.1），生成网格点用于后续遍历：

alpha1_range = linspace(0, 10, 100);
alpha2_range = linspace(0, 10, 100);
[Alpha1, Alpha2] = meshgrid(alpha1_range, alpha2_range);
bistable_map = zeros(size(Alpha1)); % 存储双稳态区域标记

2. 遍历参数，求解平衡点并判断稳定性 #

对每个(α₁, α₂)组合，求解零倾线交点（即方程组 u = α₁/(1+v^β) 和 v = α₂/(1+u^γ) 的解），再通过雅可比矩阵判断每个平衡点的稳定性：

beta = 3; gamma = 3; % 固定协同性参数
options = optimset('Display','off'); % 关闭求解过程输出

for i = 1:length(alpha1_range)
    for j = 1:length(alpha2_range)
        a1 = Alpha1(i,j);
        a2 = Alpha2(i,j);
        % 定义待求解的零倾线方程组
        fun = @(x) [x(1) - a1/(1 + x(2)^beta); x(2) - a2/(1 + x(1)^gamma)];
        % 尝试多个初始值，确保找到所有可能的平衡点（切换开关最多3个）
        initial_guesses = [[0.1;0.1], [a1;0.1], [0.1;a2]];
        eq_points = [];
        for guess = initial_guesses
            sol = fsolve(fun, guess, options);
            % 验证解的合理性（误差小于阈值，且非负）
            if norm(fun(sol)) < 1e-6 && all(sol >= 0)
                % 去重，避免重复求解相同平衡点
                if isempty(eq_points) || min(norm(eq_points - sol, 2, 2)) > 1e-3
                    eq_points = [eq_points, sol];
                end
            end
        end
        % 统计稳定平衡点数量
        stable_count = 0;
        for k = 1:size(eq_points,2)
            u_eq = eq_points(1,k);
            v_eq = eq_points(2,k);
            % 计算雅可比矩阵
            J = [ -1,  -a1*beta*v_eq^(beta-1)/(1+v_eq^beta)^2;
                  -a2*gamma*u_eq^(gamma-1)/(1+u_eq^gamma)^2,  -1 ];
            % 求解特征值判断稳定性
            eig_vals = eig(J);
            % 稳定平衡点：所有特征值实部小于0
            if all(real(eig_vals) < 0)
                stable_count = stable_count + 1;
            end
        end
        % 标记双稳态区域（存在至少2个稳定平衡点）
        if stable_count >= 2
            bistable_map(i,j) = 1;
        end
    end
end

3. 绘制分岔图 #

用填充等高线图可视化双稳态区域，同时可叠加分岔边界：

figure;
% 填充双稳态区域
contourf(Alpha1, Alpha2, bistable_map, [0.5, 1.5], 'LineColor','none');
colormap([0.8 0.8 0.8; 1 1 1]); % 灰色标记双稳态区域
xlabel('α₁'); ylabel('α₂');
title('切换开关系统双稳态区域（α₁ vs α₂）');
colorbar;
hold on;
% 绘制分岔边界
[C,h] = contour(Alpha1, Alpha2, bistable_map, [0.5, 0.5], 'Color','k','LineWidth',1.5);
clabel(C,h,'LabelSpacing',200);
hold off;

关键细节说明 #

• 初始值选择：切换开关的零倾线最多有3个交点，必须尝试多个初始值（近原点、近α₁轴、近α₂轴）才能覆盖所有可能的平衡点。
• 雅可比矩阵推导：基于系统微分方程的偏导数生成，对应：
  • ∂(du/dt)/∂u = -1，∂(du/dt)/∂v = -α₁βv^(β-1)/(1+vβ)²
  • ∂(dv/dt)/∂u = -α₂γu^(γ-1)/(1+uγ)²，∂(dv/dt)/∂v = -1
• 稳定性判断：稳定平衡点的特征值实部全为负；鞍点有一个正实部、一个负实部；双稳态区域对应存在2个稳定点+1个鞍点的参数组合。

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