如何在Matlab Simbiology中实现含时变k1的随机模拟？

问题背景 #

需要对基因表达模型进行随机模拟，其中转录速率常数k1随时间线性变化（例如k1(t) = k0 + a*t，k0为初始值，a为斜率），但现有工具默认不支持时变速率的随机模拟，已尝试Python的stochpy、gillespy2及Matlab的SimBiology工具。

解决方案：利用SimBiology事件实现参数时变更新 #

标准Gillespie算法假设速率常数恒定，要实现时变参数，可通过**事件（Event）**在模拟过程中定期更新k1的值，以此近似线性变化。具体步骤及修改后的代码如下：

% Create a new model
model = sbiomodel('GeneExpressionModel');

% Add reactions for gene expression
reaction1 = addreaction(model, 'Gene -> Gene + mRNA');
reaction2 = addreaction(model, 'mRNA -> mRNA + Protein');
reaction3 = addreaction(model, 'mRNA -> null');
reaction4 = addreaction(model, 'Protein -> null');

% Set reaction rate constants
kineticLaw1 = addkineticlaw(reaction1, 'MassAction');
kineticLaw2 = addkineticlaw(reaction2, 'MassAction');
kineticLaw3 = addkineticlaw(reaction3, 'MassAction');
kineticLaw4 = addkineticlaw(reaction4, 'MassAction');

% 将k1改为全局参数，方便事件修改
p1 = addparameter(model, 'k1', 'Value', 0.1); % 初始值k0=0.1
p2 = addparameter(kineticLaw2, 'k2', 'Value', 0.05);
p3 = addparameter(kineticLaw3, 'k3', 'Value', 0.01);
p4 = addparameter(kineticLaw4, 'k4', 'Value', 0.002);

% Set the rate equation for each reaction
kineticLaw1.ParameterVariableNames = {p1.Name};
kineticLaw2.ParameterVariableNames = {p2.Name};
kineticLaw3.ParameterVariableNames = {p3.Name};
kineticLaw4.ParameterVariableNames = {p4.Name};

% Set initial conditions
model.Species(1).InitialAmount = 1;
model.Species(2).InitialAmount = 0;
model.Species(3).InitialAmount = 0;

% --------------------------
% 添加事件实现k1的线性变化：k1(t) = 0.1 + 0.001*t
% --------------------------
% 定义斜率参数
slope_k1 = addparameter(model, 'slope_k1', 'Value', 0.001);
% 创建事件：定期触发更新k1
event = addevent(model);
% 触发规则：每隔0.1时间单位更新一次（可根据精度调整）
event.Trigger = 'time >= lastTriggerTime + 0.1';
% 更新动作：按线性公式修改k1
event.Effect = 'k1 = 0.1 + slope_k1*time';
% 设置初始触发时间
event.InitialTriggerTime = 0;

% Configure the stochastic solver
configSet = getconfigset(model, 'active');
set(configSet, 'SolverType', 'ssa'); % 启用Gillespie随机算法
set(configSet, 'StopTime', 100); % 设置模拟结束时间
% 启用事件功能
set(configSet, 'Events', 'on');

% Run the stochastic simulation
simulationData = sbiosimulate(model);

% Plot the simulation results
figure;
plot(simulationData.Time, simulationData.Data);
xlabel('Time');
ylabel('Species Count');
legend('Gene', 'mRNA', 'Protein');
title('Stochastic Simulation of Gene Expression with Time-Varying k1');

% 绘制k1的变化曲线，验证效果
figure;
t = simulationData.Time;
k1_theory = 0.1 + 0.001*t;
plot(t, k1_theory, 'b-', 'LineWidth', 1.5);
xlabel('Time');
ylabel('k1 Value');
title('Time-Varying k1 (Linear)');

关键说明 #

• 事件触发的时间间隔越小，k1的变化越接近连续线性，模拟精度越高，但会增加计算量
• 若需要严格连续的时变速率，也可直接定义反应速率为时间的函数（例如kineticLaw1.Rate = '(0.1 + 0.001*time)*Gene'），但部分SimBiology版本对SSA模式下的自定义时变速率支持有限，需自行验证
• 事件方法是SimBiology中处理时变参数的标准方案，兼容性更好

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