How do I resolve, error in resample ?

Question

0 个投票

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Follow Question

Answer 1

Rick Paxson 2022-5-30

编辑：Rick Paxson 2022-5-30

在 MATLAB Online 中打开

1 个投票

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

RAMYA BODDEPALLI 2022-6-1

Thank you so much!

请先登录，再进行评论。

Answer 2

Jeremy Huard 2022-5-30

在 MATLAB Online 中打开

1 个投票

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

RAMYA BODDEPALLI 2022-6-1

Thank you!

请先登录，再进行评论。

How do I resolve, error in resample ?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

社区

类别

产品

版本

标签

Community Treasure Hunt

How do I resolve, error in resample ?

0 个评论 显示 -2更早的评论 隐藏 -2更早的评论

采纳的回答

1 个评论 显示 -1更早的评论 隐藏 -1更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论 隐藏 -1更早的评论

社区

类别

产品

版本

标签

另请参阅

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论