Детерминированная симуляция модели, содержащей прерывистость

В этом примере показано, как правильно создать модель SimBiology ®, которая содержит разрывы.

Фон

Модель, которую вы создаете в этом примере, моделирует устранение белка первого порядка с заданной скоростью. Темп производства содержит две разрывы. Чтобы точно симулировать модель, вы должны создать события, которые запускаются во время разрыва.

Темп производства имеет три «режима» производства, как проиллюстрировано на следующем графике:

plot([0 3 3 6 6 10], [5 5 3 3 0 0]);
ylim([-0.5 5.5]);
xlabel('Time');
ylabel('Rate');
title('Discontinuous Protein Production Rate');

Figure contains an axes. The axes with title Discontinuous Protein Production Rate contains an object of type line.

Первоначально («Mode 1»), скорость производства является постоянным значением 5. От 3 до 6 секунд («Mode 2»), скорость производства 3. Через 6 секунд («Mode 3») скорость производства составляет 0. Эти скорости производства реализованы в функции MATLAB discontSimBiologyRateFunction.m, которая требует двух аргументов, времени симуляции и режима производства.

В этом примере вы добавите события в модель, чтобы изменить режим производства белка. Этот подход гарантирует, что разрывы в модели происходят только через события, что дополнительно гарантирует, что решатель ОДУ точно вычисляет численное поведение вблизи разрывов.

Обратите внимание, что для точного моделирования модели необходимо использовать события для обработки любого разрыва, независимо от того, связаны ли они со значениями функций или их производными.

Построение модели, отсека и видов

model = sbiomodel('discontinuous rate');
central = addcompartment(model,'Central');
addspecies(central,'protein')
ans = 
   SimBiology Species Array

   Index:    Compartment:    Name:      Value:    Units:
   1         Central         protein    0               

Создайте реакцию для устранения первого порядка

reaction1 = addreaction(model,'protein -> null')
reaction1 = 
   SimBiology Reaction Array

   Index:    Reaction:      
   1         protein -> null

ke = addparameter(model,'ke', 0.5);
kineticLaw1 = addkineticlaw(reaction1,'MassAction');
kineticLaw1.ParameterVariableNames = {ke.Name};
reaction1.ReactionRate;

Создайте события, которые инициируются во время разрывов

Эти события устанавливают режим параметра, который управляет режимом производства белка. Режим первоначально равен 1, изменяется на 2 во время 3 и изменяется на 3 во время 6.

counter = addparameter(model,'mode', 1, 'ConstantValue', false);
addevent(model,'time > 3', 'mode = 2')
ans = 
   SimBiology Event Array

   Index:    Trigger:    EventFcns:
   1         time > 3    mode = 2  

addevent(model,'time > 6', 'mode = 3')
ans = 
   SimBiology Event Array

   Index:    Trigger:    EventFcns:
   1         time > 6    mode = 3  

Создайте реакцию для производства белка

Мы присваиваем эту ставку параметру, используя повторное правило присвоения. Это позволяет нам хранить скорость производства в результатах симуляции.

reaction2 = addreaction(model, 'null -> protein');
rate2 = addparameter(model,'rate2', 0, 'ConstantValue', false);
reaction2.ReactionRate = 'rate2'
reaction2 = 
   SimBiology Reaction Array

   Index:    Reaction:      
   1         null -> protein

addrule(model,'rate2 = discontSimBiologyRateFunction(time, mode)', 'repeatedAssignment')
ans = 
   SimBiology Rule Array

   Index:    RuleType:             Rule:                                            
   1         repeatedAssignment    rate2 = discontSimBiologyRateFunction(time, mode)

Просмотр содержимого функции discontSimBiologyRateFunction

type discontSimBiologyRateFunction
function rate = discontSimBiologyRateFunction(time, mode)
%discontSimBiologyRateFunction - Helper function for discontSimBiologyModel demo
% RATE = discontSimBiologyRateFunction(TIME, MODE);

%   Copyright 2010 The MathWorks, Inc.

% Mode is a double precision number subject to round-off errors. We need to
% round to the nearest integer to correctly handle this issue.
mode = round(mode);
switch mode
    case 1
        rate = 5;
    case 2
        rate = 3;
    case 3
        rate = 0;
    otherwise
        error('Invalid mode.');
end
    

Симулируйте и постройте график модели

model
model = 
   SimBiology Model - discontinuous rate 

   Model Components:
     Compartments:      1
     Events:            2
     Parameters:        3
     Reactions:         2
     Rules:             1
     Species:           1
     Observables:       0

sd = sbiosimulate(model);
plot(sd.Time, sd.Data);
ylim([-0.5 8]);
xlabel('Time');
ylabel('State');
title('Simulation Results');
legend(sd.DataNames);

Figure contains an axes. The axes with title Simulation Results contains 3 objects of type line. These objects represent protein, mode, rate2.

Заключение

Этот пример иллюстрирует, как создать модель SimBiology, которая содержит разрывы. Он иллюстрирует, как добавить события в модель для устранения разрывов, чтобы можно было точно симулировать модель.