This example shows how to create a custom function and incorporate it in model simulation.
This example assumes you have a working knowledge of:
Creating and saving MATLAB programs
This example uses the model described in Model of the Yeast Heterotrimeric G Protein Cycle.
This table shows the reactions used to model the G protein cycle and the corresponding rate parameters (rate constants) for each reaction. For reversible reactions, the forward rate parameter is listed first.
|2||Heterotrimeric G protein formation|
|3||G protein activation|
|4||Receptor synthesis and degradation|
|6||G protein inactivation|
|1 Legend of species: L = ligand (alpha factor), R = alpha-factor receptor, Gd = inactive G-alpha-GDP, Gbg = free levels of G-beta:G-gamma complex, G = inactive Gbg:Gd complex, Ga = active G-alpha-GTP|
This example assumes that:
An inhibitor (
Inhib species) slows the inactivation
of the active G protein (reaction 6 above,
The variation in the amount of inhibitor (
species) is defined in a custom function,
The inhibitor (
Inhib species) affects the reaction
by changing the amount of rate parameter
This example shows how to create and call a custom function in a SimBiology® expression. Specifically, it shows how to use a custom function in a rule expression.
You can use custom functions in:
Reaction rate expressions (
Rule expressions (
Event expressions (
EventFcns property or
The requirements for using custom functions in SimBiology expressions are:
Create a custom function. For more information, see
Change the current folder to the folder containing your custom
MATLAB file. Do this by using the
cd command or by using
the Current Folder field in the MATLAB desktop toolbar. Alternatively,
add the folder containing your file to the search path. Do this by using
addpath command or see
Change Folders on the Search Path (MATLAB).
Call the custom function in a SimBiology reaction, rule, or event expression.
If your rule or reaction rate expression is not continuous and differentiable, see Using Events to Address Discontinuities in Rule and Reaction Rate Expressions before simulating your model.
The following procedure creates a custom function,
inhibvalex, which lets you specify how the inhibitor amount
changes over time. The inputs are time, the initial amount of inhibitor, and a
parameter that governs the amount of inhibitor. The output of the function is the
amount of inhibitor.
In the MATLAB desktop, select File > New > Script, to open the MATLAB Editor.
Copy and paste the following function declaration:
% inhibvalex.m function Cp = inhibvalex(t, Cpo, kel) % This function takes the input arguments t, Cpo, and kel % and returns the value of the inhibitor Cp. % You can later specify the input arguments in a % SimBiology rule expression. % For example in the rule expression, specify: % t as time (a keyword recognized as simulation time), % Cpo as a parameter that represents the initial amount of inhibitor, % and kel as a parameter that governs the amount of inhibitor. if t < 400 Cp = Cpo*exp(-kel*(t)); else Cp = Cpo*exp(-kel*(t-400)); end
Save the file (name the file
inhibvalex.m) in a
directory that is on the MATLAB search path, or to a directory that you can access.
If the location of the file is not on the MATLAB search path, change the working directory to the file location.
gprotein example project, which includes the variable
m1, a model object:
m1 model object appears in the MATLAB Workspace.
The following procedure creates a rule expression that calls the custom function,
inhibvalex, and specifies the three input values to this
Add a repeated assignment rule to the model that specifies the three input
values to the custom function,
rule1 = addrule(m1, 'Inhib = inhibvalex(time, Cpo, Kel)',... 'repeatedAssignment');
time input is a SimBiology keyword recognized as simulation time
Create the two parameters used by the
rule1 rule and
assign values to them:
p1 = addparameter(m1, 'Cpo', 250); p2 = addparameter(m1, 'Kel', 0.01);
Create the species used by the
s1 = addspecies(m1.Compartments, 'Inhib');
The value of rate parameter
kGd is affected by the amount of
inhibitor present in the system. Add a rule to the model to describe this action,
but first change the
ConstantValue property of the parameter
kGd so that it can be varied by a rule.
ConstantValue property of the
kGd parameter to
p3 = sbioselect(m1, 'Type', 'parameter', 'Name', 'kGd'); p3.ConstantValue = false;
Add a repeated assignment rule to the model to define how the
kGd parameter is affected by the
rule2 = addrule(m1, 'kGd = 1/Inhib', 'repeatedAssignment');
The custom function,
inhibvalex, introduces a discontinuity
in the model when time = 400. To ensure accurate simulation results, add an event to
the model to reset the solver at the time of the discontinuity. Set the event to
trigger at the time of the discontinuity (time = 400). The event does not need to
modify the model, so create an event function that multiplies a species value by
addevent(m1, 'time>=400', 'G=1*G');
Configure the simulation settings (
m1 model object to log all states during the
cs = getconfigset(m1); cs.RuntimeOptions.StatesToLog = 'all';
Simulate the model.
simDataObj = sbiosimulate(m1);
Plot the results.
The plot does not show the species of interest due to the wide range in species amounts/concentrations.
Plot only the species of interest.
GaSimDataObj = selectbyname(simDataObj,'Ga'); sbioplot(GaSimDataObj);
Notice the change in the profile of species
Ga at time
400 seconds (simulation time). This is the time when
the inhibitor amount is changed to reflect the re-addition of inhibitor to
Plot only the inhibitor (
InhibSimDataObj = selectbyname(simDataObj,'Inhib'); sbioplot(InhibSimDataObj)