createSimFunction does not use initial-as​signment-c​omputed values for MassAction ParameterVariableNames

When a MassAction rate parameter (e.g., `ke`) is set by an initial assignment rule, `createSimFunction` does not use the computed value for the MassAction rate calculation. The parameter IS correctly computed (confirmed by observing it as a SimFunction output), but MassAction uses the pre-assignment `.Value` property instead.

`sbiosimulate` with the same model produces correct results.

Environment

- MATLAB R2025b (25.2.0.3123386)

- SimBiology toolbox

model = sbiomodel("MWE");

cs = getconfigset(model);
cs.SolverType = 'sundials';
cs.SolverOptions.AbsoluteTolerance = 1e-8;
cs.SolverOptions.RelativeTolerance = 1e-8;
cs.CompileOptions.DimensionalAnalysis = true;
cs.CompileOptions.UnitConversion = true;

addcompartment(model, 'Central', ...
    'CapacityUnits', 'milliliter', 'ConstantCapacity', false);

addspecies(model.Compartments(1), 'Drug', 0, 'Units', 'milligram/milliliter');

addparameter(model, 'typical_V', 3000, 'Units', 'milliliter');
addparameter(model, 'typical_CL', 200, 'Units', 'milliliter/hour');
addparameter(model, 'BW', 70, 'Units', 'kilogram');
addparameter(model, 'typical_BW', 70, 'Units', 'kilogram');
addparameter(model, 'CL', eps, 'Units', 'milliliter/hour');
addparameter(model, 'ke', eps, 'Units', '1/hour');

% Initial assignments:
%   Central = 3000 mL, CL = 200 mL/h, ke = 200/3000 = 0.0667 /h
addrule(model, 'Central = typical_V * (BW / typical_BW)', 'initialAssignment');
addrule(model, 'CL = typical_CL * (BW / typical_BW)', 'initialAssignment');
addrule(model, 'ke = CL / Central', 'initialAssignment');

% MassAction elimination (first-order)
rxn = addreaction(model, 'Central.Drug -> null', 'Name', 'Elimination');
kl = addkineticlaw(rxn, 'MassAction');
kl.ParameterVariableNames = 'ke';

% Dose
dose = sbiodose('IV', 'schedule');
dose.TargetName = 'Drug';
dose.Amount = 300;
dose.AmountUnits = 'milligram';
dose.TimeUnits = 'hour';
dose.Time = 0;

% --- Method 1: createSimFunction (INCORRECT) ---
F = createSimFunction(model, {}, {'Drug', 'ke', 'CL', 'Central'}, {'Central.Drug'}, ...
    'UseParallel', false, 'AutoAccelerate', false);
[~, y] = F([], [], getTable(dose), [0, 10, 100]);

fprintf('createSimFunction:\n');
fprintf('  Drug(t=0)=%g, Drug(t=10h)=%g, Drug(t=100h)=%g\n', y{1}(1,1), y{1}(2,1), y{1}(3,1));
fprintf('  ke=%g, CL=%g, Central=%g  (all correct!)\n\n', y{1}(1,2), y{1}(1,3), y{1}(1,4));

% --- Method 2: sbiosimulate (CORRECT) ---
cs.StopTime = 100;
cs.TimeUnits = 'hour';
cs.SolverOptions.OutputTimes = [0, 10, 100];
[~, x, names] = sbiosimulate(model, cs, [], dose);
idx = strcmp(names, "Drug");

fprintf('sbiosimulate:\n');
fprintf('  Drug(t=0)=%g, Drug(t=10h)=%g, Drug(t=100h)=%g\n', x(1,idx), x(2,idx), x(3,idx));

Output

createSimFunction:
Drug(t=0)=0.1, Drug(t=10h)=0.0999815, Drug(t=100h)=0.099815
ke=0.0666667, CL=200, Central=3000  (all correct!)

sbiosimulate:
Drug(t=0)=0.1, Drug(t=10h)=0.0513417, Drug(t=100h)=0.000127269

Expected Behavior

Both methods should produce identical results. Expected half-life = ln(2) / 0.0667 = 10.4 hours. `sbiosimulate` is correct; `createSimFunction` shows ~1000x slower elimination.

Key Observation

The initial assignments ARE evaluated correctly in `createSimFunction` — `ke=0.0667`, `CL=200`, `Central=3000` are all reported with correct values when observed as SimFunction outputs. However, the MassAction rate calculation does not use the computed `ke` value; it appears to use the pre-assignment `.Value` (eps).

Notes

- This pattern (`ke = CL/V` as initial assignment, MassAction with `ke`) is the same pattern used by SimBiology's own `PKCompartment.m` source code for "linear-clearance" elimination (lines 132-142 in R2025b).

- The issue is not specific to chained initial assignments — even `ke = typical_CL / typical_V` (single-level, depending only on directly-set parameters) exhibits the same behavior.

- The `Constant` property of the parameters does not affect the outcome.

- `AutoAccelerate = true` does not fix the issue.

Answer 1

[Jeremy Huard]

Hi @Paulo,

the discrepancy comes from the time units used in both cases.

You simulate the SimFunction with Time units in second but then use time units of hour for sbiosimulate.

If you add the line

cs.TimeUnits = 'hour';

before creating the SimFunction, you will get the same results in both cases.

Best,

Jérémy

[Paulo]

Hi Jeremy, you are right, that solves the problem.

I imagined that the DimensionalAnalysis and UnitConversion both set to True would manage this situation.

Thanks for the help !

[Jeremy Huard]

The SimFunction F you created has a TimeUnits property, which corresponds to the configset TimeUnits property at the time of creation.

So, when you provide output times or a stop time to simulate the SimFunction, make sure the numerical values are given in that unit.