simbio.complexstep.max

max function for SimBiology local sensitivity analysis

Since R2022b

Syntax

M = simbio.complexstep.max(X,Y)

Description

M = simbio.complexstep.max(X,Y) returns X if real(X) is greater than real(Y). Otherwise, the function returns Y. simbio.complexstep.max is equivalent to max when both inputs are real.

• When you have a SimBiology model with a custom function that calls max and you are performing local sensitivity analysis on the model, replace max with simbio.complexstep.max in your custom function.

• You do not need to update SimBiology expressions (such as reaction rates or rules) that directly call max. SimBiology automatically replaces max with simbio.complexstep.max whenever:

  • You calculate local sensitivities using sbiosimulate, a SimFunctionSensitivity object, or the Calculate Sensitivities program.

  • sbiofit or fitproblem uses local sensitivity analysis to determine the gradients of the objective function during parameter estimation.

• simbio.complexstep.max is not differentiable when the real parts of X and Y are equal [1]. For more information, see Local Sensitivity Analysis (LSA).

Examples

Use Replacements for abs, min, and max in Custom Functions for SimBiology Local Sensitivity Analysis

This example shows how to use replacements for abs, min, and max in some custom functions so that the model becomes compatible with local sensitivity analysis (LSA). The replacement functions are simbio.complexstep.abs, simbio.complexstep.min, and simbio.complexstep.max. Specifically, one of the custom functions used by this example computes the net amount of a drug species moved between two compartments. The other custom function sets the thresholds for the forward and reverse reaction fluxes.

Create a model.

model       = sbiomodel("model");
c1          = addcompartment(model,"c1");
c2          = addcompartment(model,"c2");
s1          = addspecies(c1,"Drug",2);
s2          = addspecies(c2,"Drug");
netFlux     = addspecies(c1,"netFlux");
reaction    = addreaction(model,"c1.Drug <-> c2.Drug");
kf          = addparameter(model,"kf",1.0);
kr          = addparameter(model,"kr",1.5);
fluxMin     = addparameter(model,"fluxMin",0.1);
fluxMax     = addparameter(model,"fluxMax",10);

Define the net amount of drug species moved between two compartments using a custom function calculateNetFlux based on the constrained forward and reverse reaction fluxes, which are defined later.

rule = addrule(model,"netFlux = calculateNetFlux(boundedForwardFlux,boundedReverseFlux)","rate");

calculateNetFlux uses simbio.complexstep.abs, and the function is already saved in the provided file named calculateNetFlux.m.

type calculateNetFlux.m

function netFlux = calculateNetAmount(forwardFlux,reverseFlux)
    netFlux   = simbio.complexstep.abs(forwardFlux-reverseFlux);
end

Define the forward and reverse fluxes of the reaction. Set the thresholds on the fluxes using the imposeBounds custom function.

boundedForwardFlux = addparameter(model,"boundedForwardFlux","Constant",false);
boundedReverseFlux = addparameter(model,"boundedReverseFlux","Constant",false);
forwardFlux        = addparameter(model,"forwardFlux","Constant",false);
reverseFlux        = addparameter(model,"reverseFlux","Constant",false);

forwardFlux        = addrule(model,"forwardFlux = kf*c1.Drug","repeatedAssignment");
reverseFlux        = addrule(model,"reverseFlux = kr*c2.Drug","repeatedAssignment");
boundedForwardFlux = addrule(model,"boundedForwardFlux = imposeBounds(forwardFlux,fluxMin,fluxMax)","repeatedAssignment");
boundedReverseFlux = addrule(model,"boundedReverseFlux = imposeBounds(reverseFlux,fluxMin,fluxMax)","repeatedAssignment");
reaction.ReactionRate = "boundedForwardFlux - boundedReverseFlux";

imposeBounds uses simbio.complexstep.min and simbio.complexstep.max to set the lower and upper limits for the reaction flux.

type imposeBounds.m

function boundedFlux = imposeBounds(fluxInput,fluxMin,fluxMax)
    fm             = simbio.complexstep.max(fluxMin,fluxInput);
    boundedFlux    = simbio.complexstep.min(fluxMax,fm);
end

Enable local sensitivity analysis.

configset = getconfigset(model);
configset.RuntimeOptions.StatesToLog        = "netFlux";
configset.SolverOptions.SensitivityAnalysis = true;
sensitivityOptions          = configset.SensitivityAnalysisOptions;
sensitivityOptions.Inputs   = [kf,kr];
sensitivityOptions.Outputs  = netFlux;

Temporarily disable the warning about unsupported functions. The warning is safe to ignore.

warnState  = warning("off","SimBiology:senscsverify:UnsupportedFunction");
cleanupobj = onCleanup(@()warning(warnState));

Simulate the model and perform sensitivity analysis.

simdata = sbiosimulate(model);
sbioplot(simdata);

delete(cleanupobj);

Input Arguments

X — First input

scalar

First input, specified as a numeric scalar.

Data Types: double
Complex Number Support: Yes

Y — Second input

scalar

Second input, specified as a numeric scalar.

Data Types: double
Complex Number Support: Yes

References

[1] Martins, Joaquim, Ilan Kroo, and Juan Alonso. “An Automated Method for Sensitivity Analysis Using Complex Variables.” In 38th Aerospace Sciences Meeting and Exhibit. Reno, NV, U.S.A.: American Institute of Aeronautics and Astronautics, 2000. https://doi.org/10.2514/6.2000-689.

Version History

Introduced in R2022b