Using 'time' in observable expressions

Hello,

Is it possible to select only certain timepoints for observable calculations, for example:

mdl.addobservable('maxA','max(A(time>(3*7*24)))','Units','dimensionless');

In my example, A is a parameter that is set by a repeated assignment rule.

When I run the model in R2021b, I get the following error, whih suggests incorrect dimensions (but my observable expression should return a scalar, no?):

------------------------------------------------------

When running sbiosteadystate:

Error using sbiosteadystate (line 128)

sbiosteadystate encountered a model compilation error.

Caused by:

Error using SimBiology.internal.verifyHelper

--> Error reported from Expression Validation:

The result of evaluating observable 'maxA' has the wrong size. It must either be a scalar or a column vector of the same length as the time vector.

-------------------------------------------------------

When running sbiosimulate:

Error using SimBiology.internal.simulate

--> Error reported from Expression Validation:

The result of evaluating observable 'maxA' has the wrong size. It must either be a scalar or a column vector of the same length as the time vector.

----------------------------------------------------

Please advise.

Thank you,

Abed

Answer 1

[Florian Augustin]

Hi Abed,

when evaluating expressions of observables, SimBiology first tries to determine whether the result is a scalar value or a vector. I therefore recommend to write observables' expressions in a way that they either always return a scalar value or always a vector of values - for any possible input values of referenced model components and time.

In your case, the expression max(A(time > 3*7*24))) can be a vector whose length differs from the length of the time vector, e.g. if none of the values in time exceeds 3*7*24, or if the maximum is attained at multiple time points. You can write a helper function that returns the max. value and protects against invalid edge cases:

function maxValue = getMaxValueAfterSpecifiedTime(values, time, threshold)

% Get max value (call unique to ensure scalar or empty output).

maxValue = unique(max(values(time > threshold)));

% Protect against empty maxValue by returning nan.

if isempty(maxValue)

maxValue = nan;

end

end

You can then call this helper as follows:

% Load/build your model

mdl = sbmlimport("radiodecay");

% Observable referencing time point

observable = addobservable(mdl, "maxX", "getMaxValueAfterSpecifiedTime(x, time, 2)");

set(observable, "Active", true);

% Set simulation options

configset = getconfigset(mdl);

configset.RuntimeOptions.StatesToLog = "x";

% Simulate and plot model state and active observable

simData = sbiosimulate(mdl);

sbioplot(simData);

When you run this particular example, you may see that maxX differs from the actual max value as shown in the plots. You can add interpolation to your helper function to increase the accuracy of the reported max value; an example how this could look like is attached to this post.

I hope this helps.

Best,

-Florian

[Abed Alnaif]

Thank you, Florian! This makes sense and this was very helpful. Out of curiosity, is it required to set the Active property of the observable to true?

[Florian Augustin]

Hi Abed,

I'm glad it helped!

The line where I set the Active property on the observable is actually not necessary, the default value for this property is already true. I only included it to make it obvious at first sight what to expect to see in the plot when running the code (active observables and states (values of species, parameters, compartments) listed in the configuration set's StatesToLog property).

Best,

Florian

[Arthur Goldsipe]

I just had an idea for how to avoid writing a separate function getMaxValueAfterSpecifiedTime. In many cases, I think you could instead change

max(expression)

to

max(vertcat(nan,expression))

to make sure that max is never called on an empty vector.

[Abed Alnaif]

Thanks for sharing, Arthur. I've been using the getMaxValueAfterSpecifiedTime function and it's been working, but this is is a neat trick.