How do I resolve, error in resample ?

8 views (last 30 days)
I am trying to model a positive autoregulation network.
I ran it for 1000 ensemble runs, while trying to get the statistics, it shows the following error. Could you please help me with this?
model = sbiomodel('PAR Model');
c = addcompartment(model,'C');
c.CapacityUnits = 'meter^3';
Enter Reactions
r1 = addreaction(model,'off + P -> on + P')
r1 =
SimBiology Reaction Array Index: Reaction: 1 off + P -> on + P
r2 = addreaction(model,'on -> off')
r2 =
SimBiology Reaction Array Index: Reaction: 1 on -> off
r3 = addreaction(model,'on -> on + R ')
r3 =
SimBiology Reaction Array Index: Reaction: 1 on -> on + R
r4 = addreaction(model, 'R -> null')
r4 =
SimBiology Reaction Array Index: Reaction: 1 R -> null
r5 = addreaction(model, 'R -> R + P')
r5 =
SimBiology Reaction Array Index: Reaction: 1 R -> R + P
r6 = addreaction(model, 'P -> null')
r6 =
SimBiology Reaction Array Index: Reaction: 1 P -> null
Set Reactions to be MassAction
k1 = addkineticlaw(r1, 'MassAction');
k2 = addkineticlaw(r2, 'MassAction');
k3 = addkineticlaw(r3, 'MassAction');
k4 = addkineticlaw(r4, 'MassAction');
k5 = addkineticlaw(r5, 'MassAction');
k6 = addkineticlaw(r6, 'MassAction');
Add Rate Constants for Each Reaction
p1 = addparameter(k1, 'kpar', 'Value', 1 , 'ValueUnits', '1/(minute*molecule)');
p2 = addparameter(k2, 'koff', 'Value', 28, 'ValueUnits', '1/minute' );
p3 = addparameter(k3, 'ktransc', 'Value', 1.2, 'ValueUnits', '1/minute' );
p4 = addparameter(k4, 'kdegR', 'Value', 0.36,'ValueUnits', '1/minute' );
p5 = addparameter(k5, 'ktransl', 'Value', 6,'ValueUnits', '1/minute' );
p6 = addparameter(k6, 'kdegP', 'Value', 2.5, 'ValueUnits', '1/minute' );
Set the Kinetic Law Constants for Each Kinetic Law.
k1.ParameterVariableNames = {'kpar'};
k2.ParameterVariableNames = {'koff'};
k3.ParameterVariableNames = {'ktransc'};
k4.ParameterVariableNames = {'kdegR'};
k5.ParameterVariableNames = {'ktransl'};
k6.ParameterVariableNames = {'kdegP'};
Specify Initial Amounts of Each Species
model.species(1).InitialAmount = 1 ;
model.species(1).InitialAmountUnits = 'molecule' % off
model =
SimBiology Model - PAR Model Model Components: Compartments: 1 Events: 0 Parameters: 6 Reactions: 6 Rules: 0 Species: 4 Observables: 0
model.species(2).InitialAmount = 25 ;
model.species(2).InitialAmountUnits = 'molecule' % P
model =
SimBiology Model - PAR Model Model Components: Compartments: 1 Events: 0 Parameters: 6 Reactions: 6 Rules: 0 Species: 4 Observables: 0
model.species(3).InitialAmount = 0 ;
model.species(3).InitialAmountUnits = 'molecule' % on
model =
SimBiology Model - PAR Model Model Components: Compartments: 1 Events: 0 Parameters: 6 Reactions: 6 Rules: 0 Species: 4 Observables: 0
model.species(4).InitialAmount = 1 ;
model.species(4).InitialAmountUnits = 'molecule' % R
model =
SimBiology Model - PAR Model Model Components: Compartments: 1 Events: 0 Parameters: 6 Reactions: 6 Rules: 0 Species: 4 Observables: 0
Get the Active Configuration Set for the Model.
cs = getconfigset(model,'active');
Simulate Model Using SSA Stochastic Solver
numRuns = 1000;
cs.SolverType = 'ssa';
cs.SolverOptions.LogDecimation = 100;
simdata = sbioensemblerun(model, numRuns);
figure;
speciesNames = {'R'};
[t, x] = selectbyname(simdata, speciesNames);
hold on;
for i = 1:numRuns
plot(t{i}, x{i});
end
grid on;
xlabel('Time in seconds');
ylabel('Species amount');
title('Raw Ensemble Data');
legend(speciesNames);
[timeVals, meanVals, varianceVals] = sbioensemblestats(simdata, speciesNames);
Error using SimData/resample
SimData objects must have at least two time points to be resampled.

Error in sbioensemblestats (line 100)
filtersd = filtersd.resample([],interp);

Accepted Answer

Ricardo Paxson
Ricardo Paxson on 30 May 2022
Edited: Ricardo Paxson on 30 May 2022
Thank you for reporting this. It appears to be a bug in sbioensemblestats. The default interpolation method used requires at least 2 points for each simulation profile and your example does not satisfly that requirement. You can see this by looking at the Time property of the simdata array (e.g. >> simdata.Time). Until we fix this with a better warning/error message, I can suggest a workaround.
The decimation of the resulting simulation data is set to 100 in your example. If you set that to a lower number (I tried 10) you will get more points and the interpolation onto a common time vector will work.
cs.SolverOptions.LogDecimation = 10;
simdata = sbioensemblerun(model, numRuns);
sbioensembleplot(simdata, 'R');
Hope this helps, and thank you again for reporting this problem.

More Answers (1)

Jeremy Huard
Jeremy Huard on 30 May 2022
Hi Ramya,
many of the simulations generated with this code contains only one data point. sbioensemblestats tries first to synchronize all simulation to the same time vector and it fails in this case.
In your example, this is due to the large value of LogDecimation. You could set to 1 here.
Also, UnitConversion is set to false, so your model is simulated with the current time unit in your model (minute). Which is probably what you want. But I still recommend to turn on UnitConversion and adapt the simulation stopTime accordingly:
cs.StopTime = 10;
cs.TimeUnits = 'minute';
cs.CompileOptions.UnitConversion = true;
cs.SolverType = 'ssa';
cs.SolverOptions.LogDecimation = 1;
With these settings, your code should run:
model = sbiomodel('PAR Model');
c = addcompartment(model,'C');
c.CapacityUnits = 'meter^3';
r1 = addreaction(model,'off + P -> on + P');
r2 = addreaction(model,'on -> off');
r3 = addreaction(model,'on -> on + R ');
r4 = addreaction(model, 'R -> null');
r5 = addreaction(model, 'R -> R + P');
r6 = addreaction(model, 'P -> null');
k1 = addkineticlaw(r1, 'MassAction');
k2 = addkineticlaw(r2, 'MassAction');
k3 = addkineticlaw(r3, 'MassAction');
k4 = addkineticlaw(r4, 'MassAction');
k5 = addkineticlaw(r5, 'MassAction');
k6 = addkineticlaw(r6, 'MassAction');
p1 = addparameter(k1, 'kpar', 'Value', 1 , 'ValueUnits', '1/(minute*molecule)');
p2 = addparameter(k2, 'koff', 'Value', 28, 'ValueUnits', '1/minute' );
p3 = addparameter(k3, 'ktransc', 'Value', 1.2, 'ValueUnits', '1/minute' );
p4 = addparameter(k4, 'kdegR', 'Value', 0.36,'ValueUnits', '1/minute' );
p5 = addparameter(k5, 'ktransl', 'Value', 6,'ValueUnits', '1/minute' );
p6 = addparameter(k6, 'kdegP', 'Value', 2.5, 'ValueUnits', '1/minute' );
k1.ParameterVariableNames = {'kpar'};
k2.ParameterVariableNames = {'koff'};
k3.ParameterVariableNames = {'ktransc'};
k4.ParameterVariableNames = {'kdegR'};
k5.ParameterVariableNames = {'ktransl'};
k6.ParameterVariableNames = {'kdegP'};
model.species(1).InitialAmount = 1 ;
model.species(1).InitialAmountUnits = 'molecule'; % off
model.species(2).InitialAmount = 25 ;
model.species(2).InitialAmountUnits = 'molecule'; % P
model.species(3).InitialAmount = 0 ;
model.species(3).InitialAmountUnits = 'molecule'; % on
model.species(4).InitialAmount = 1 ;
model.species(4).InitialAmountUnits = 'molecule'; % R
cs = getconfigset(model,'active');
numRuns = 1000;
cs.StopTime = 10;
cs.TimeUnits = 'minute';
cs.CompileOptions.UnitConversion = true;
cs.SolverType = 'ssa';
cs.SolverOptions.LogDecimation = 1;
simdata = sbioensemblerun(model, numRuns, cs);
figure;
speciesNames = {'R'};
[t, x] = selectbyname(simdata, speciesNames);
hold on;
for i = 1:numRuns
plot(t{i}, x{i});
end
set(gca, 'XLimitMethod','padded','YLimitMethod','padded')
grid on;
xlabel("Time in " + cs.TimeUnits);
ylabel('Species amount');
title('Raw Ensemble Data');
[timeVals, meanVals, varianceVals] = sbioensemblestats(simdata, speciesNames);
figure;
plot(timeVals, meanVals, '-', ...
timeVals, sqrt(varianceVals), 'r:', LineWidth=1.5);
grid on;
box off;
set(gca, 'XLimitMethod','padded','YLimitMethod','padded')
xlabel("Time in " + cs.TimeUnits);
title('Mean and Standard Deviation');
legend('Mean', 'Standard Deviation',Box='off');
Best regards,
Jérémy

Communities

More Answers in the  SimBiology Community

Categories

Find more on Scan Parameter Ranges in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by