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

В этом примере показано, как правильно создать модель 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.

Первоначально ("Режим 1"), производительность является постоянным значением 5. С 3 до 6 секунд ("Режим 2"), производительность равняется 3. После 6 секунд ("Режим 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, которая содержит разрывы. Это иллюстрирует, как добавить события в модель, чтобы обратиться к разрывам, таким образом, можно симулировать модель точно.