www.gusucode.com > simbio 案例源码程序 matlab代码 > simbio/discontSimBiologyModel.m

    %% Deterministic Simulation of a Model Containing a Discontinuity
%
% This example shows how to correctly build a SimBiology(R) model
% that contains discontinuities.

% Copyright 2010-2014 The MathWorks, Inc.

%% Background
% The model you create in this example simulates the first-order elimination
% of a protein that is produced at a specified rate. The production rate
% contains two discontinuities. To simulate the model accurately, you must
% create events that are triggered at the time of the discontinuity.
%
% The production rate has three "modes" of production, as illustrated in
% the following plot:
%
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');
%%
% Initially ("Mode 1"), the production rate is a constant value of 5. From
% 3 to 6 seconds ("Mode 2"), the production rate is 3. After 6 seconds
% ("Mode 3"), the production rate is 0. These production rates are
% implemented in a MATLAB function discontSimBiologyRateFunction.m, which
% requires two arguments, simulation time and production mode.
%
% In this example, you will add events to the model to change the mode of
% protein production. This approach ensures that discontinuities in the
% model occur only via events, which further ensures that the ODE solver
% accurately calculates the numerical behavior near the discontinuities.
%
% Note that to simulate a model accurately you must use events to handle
% _any_ discontinuity, whether related to function values or their
% derivatives.

%% Construct the Model, Compartment, and Species
model = sbiomodel('discontinuous rate');
central = addcompartment(model,'Central');
addspecies(central,'protein')

%% Construct the Reaction for First-Order Elimination
reaction1 = addreaction(model,'protein -> null')
ke = addparameter(model,'ke', 0.5);
kineticLaw1 = addkineticlaw(reaction1,'MassAction');
kineticLaw1.ParameterVariableNames = {ke.Name};
reaction1.ReactionRate;

%% Construct the Events That Are Triggered at the Time of Discontinuities
% These events set a parameter mode that controls the mode of protein
% production. The mode is initially 1, changes to 2 at time 3, and changes
% to 3 at time 6.
counter = addparameter(model,'mode', 1, 'ConstantValue', false);
addevent(model,'time > 3', 'mode = 2')
addevent(model,'time > 6', 'mode = 3')

%% Construct the Reaction for Protein Production
% We assign this rate to a parameter using a repeated assignment rule. This
% lets us store the production rate in the simulation results.
reaction2 = addreaction(model, 'null -> protein');
rate2 = addparameter(model,'rate2', 0, 'ConstantValue', false);
reaction2.ReactionRate = 'rate2'
addrule(model,'rate2 = discontSimBiologyRateFunction(time, mode)', 'repeatedAssignment')

%% View the Contents of discontSimBiologyRateFunction
type discontSimBiologyRateFunction

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

%% Conclusion
% This example illustrates how to create a SimBiology model that contains
% discontinuities. It illustrates how to add events to the model to address
% the discontinuities, so you can simulate the model accurately.