This example shows how to correctly build a SimBiology® model that contains discontinuities.

**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')

SimBiology Species Array Index: Compartment: Name: InitialAmount: InitialAmountUnits: 1 Central protein 0

**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;

SimBiology Reaction Array Index: Reaction: 1 protein -> null

**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')

SimBiology Event Array Index: Trigger: EventFcns: 1 time > 3 mode = 2 SimBiology Event Array Index: Trigger: EventFcns: 1 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')

SimBiology Reaction Array Index: Reaction: 1 null -> protein SimBiology Rule Array Index: RuleType: Rule: 1 repeatedAssignment rate2 = discontSimBiologyRateFunction(time, mode)

**View the Contents of 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

**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);

SimBiology Model - discontinuous rate Model Components: Compartments: 1 Events: 2 Parameters: 3 Reactions: 2 Rules: 1 Species: 1

**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.

Was this topic helpful?