How do I resolve, error in resample ?

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);

Answer 1

[Rick Paxson]

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.

[RAMYA BODDEPALLI]

Thank you so much!

Answer 2

[Jeremy Huard]

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

[RAMYA BODDEPALLI]

Thank you!