Hi Sam,
I'm sorry you're having such a hard time performing the calculations you want. Let me try to clarify things further.
I think the thing that Pramod and I are having trouble communicating to you is that you cannot just pass a vector or matrix of numbers to "sbiosimulate". You must set the parameter values some other way. One way is to change the Value property of a parameter. Here's a bit of sample code based on your code that sets the value of 'a1' to 0:
parameterObject = addparameter(kineticlaw1, 'a1',23);
parameterObject.Value = 0;
However, Pramod and I both originally thought that the best way for you to do what you want is to use something called a SimBiology variant object. A variant is a MATLAB object, just like a SimBiology reaction is an object. You use a variant to temporarily override parameter values (or species initial amounts). If you need more details, please search the SimBiology documentation. But in brief, you can create a variant using the "sbiovariant" command, which is what Pramod and I did in our examples. A variant must have a name, which is the (first) argument we pass to "sbiovariant." But to make a variant useful, you have to set its Content property. This property is a cell array that specifies the names of the parameters you want to change and their new values. You can then include the variant as argument to "sbiosimulate". Here's an example based on your code:
variant1 = sbiovariant('variant 1');
% Replace "Reaction1" with the name of your reaction.
% If the reaction has no name, replace the 1 in Reaction1 with the reaction's index.
variant1.addcontent({'parameter', 'Reaction1.a1', 'Value', 0});
sd1 = sbiosimulate(car, variant1); % Simulate model "car" with "a1=0"
You can also pass a vector of variant objects to "sbiosimulate". Pramod was suggesting that you use this capability as an easy way to loop through all the knockdown simulations you want to perform. However, because it can be awkward to construct variants for parameters that are added to kinetic laws, I'm going to show you how you can perform your knockdown simulations by setting the Value properties of the parameters, assuming your model is stored in variable "car":
allParameters = sbioselect(car, 'Type', 'parameter');
numParameters = numel(allParameters);
maxNumKnockdowns = 3;
simresultVector = SimData.empty; % Initialize the result vector
for numKnockdowns = 1:maxNumKnockdowns
allKnockdowns = nchoosek(1:numParameters, numKnockdowns);
for thisKnockdownIndex = 1:size(allKnockdowns, 1)
thisKnockdown = allKnockdowns(thisKnockdownIndex, :);
% Save the original parameter values and set them to 0.
for thisParameterIndex = 1:size(allKnockdowns, 2)
thisParameter = allParameters(thisKnockdown(thisParameterIndex));
originalParameterValues(thisParameterIndex) = thisParameter.Value;
thisParameter.Value = 0;
end
simresultVector(end+1) = sbiosimulate(car);
% Restore the original parameter values.
for thisParameterIndex = 1:size(allKnockdowns, 2)
thisParameter = allParameters(thisKnockdown(thisParameterIndex));
thisParameter.Value = originalParameterValues(thisParameterIndex);
end
end
end
You can now call "getsensmatrix" on "simresultVector(i)" to see the sensitivities for the i'th knockdown experiment.
I hope this clears things up for you.
Arthur (& Pramod)
"Sam " <sameiwater@gmail.com> wrote in message <jioqp4$e65$1@newscl01ah.mathworks.com>...
> Hello,
>
> I spent my whole day to figure out how I could insert your command lines in this script, but, I didn't get any output. Here is the latest draft of my script that may let you better appreciate where I got stuck.
>
> Paramvect=[23 13 15 0.1 143 12 50 27 34];
> %1
> kineticlaw1 = addkineticlaw(react1, 'MassAction');
> addparameter(kineticlaw1, 'a1',23);
> addparameter(kineticlaw1, 'a2',13);
> set(kineticlaw1,'ParameterVariableNames',{'a1', 'a2'});
> set(car.Reactions(1).Reactants(1), 'InitialAmount', 6);
> set(car.Reactions(1).Reactants(2), 'InitialAmount', 12);
> set(car.Reactions(1).Products(1), 'InitialAmount', 0);
> %2
> kineticlaw2 = addkineticlaw(react2,'MassAction');
> addparameter(kineticlaw2, 'a3',15);
> set(kineticlaw2,'ParameterVariableNames',{'a3'});
> set(car.Reactions(2).Reactants(1), 'InitialAmount', 0);
> set(car.Reactions(2).Products(1), 'InitialAmount', 54);
> set(car.Reactions(2).Products(2), 'InitialAmount', 12);
> %3
> kineticlaw3 = addkineticlaw(react3, 'MassAction');
> addparameter(kineticlaw3, 'a4',0.1);
> addparameter(kineticlaw3, 'a5',143);
> set(kineticlaw3,'ParameterVariableNames',{'a4', 'a5'});
> set(car.Reactions(3).Reactants(1), 'InitialAmount', 222);
> set(car.Reactions(3).Reactants(2), 'InitialAmount', 11);
> set(car.Reactions(3).Products(1), 'InitialAmount', 0);
> %4
> kineticlaw4 = addkineticlaw(react4, 'MassAction');
> addparameter(kineticlaw4, 'a6',12);
> set(kineticlaw4,'ParameterVariableNames',{'a6'});
> set(car.Reactions(4).Reactants(1), 'InitialAmount', 431);
> set(car.Reactions(4).Products(1), 'InitialAmount', 65);
> set(car.Reactions(4).Products(2), 'InitialAmount', 0);
> % 5
> kineticlaw5 = addkineticlaw(react5, 'MassAction');
> addparameter(kineticlaw5,'a7',50);
> addparameter(kineticlaw5,'a8',27);
> set(kineticlaw5,'ParameterVariableNames',{'a7', 'a8'});
> set(car.Reactions(5).Reactants(1), 'InitialAmount', 65);
> set(car.Reactions(5).Reactants(2), 'InitialAmount', 0);
> set(car.Reactions(5).Products(1), 'InitialAmount', 653);
> % 6
> kineticlaw6 = addkineticlaw(react6, 'MassAction');
> addparameter(kineticlaw6, 'a9',34);
> set(kineticlaw6,'ParameterVariableNames',{'a9'});
> set(car.Reactions(6).Reactants(1), 'InitialAmount', 0.3);
> set(car.Reactions(6).Reactants(2), 'InitialAmount', 8);
> set(car.Reactions(6).Products(1), 'InitialAmount', 0.9);
>
> I made "car" model using "Paramvect". The next step is to make a "PCar" model using a matrix of perturbedparameters called "Variant"
>
> % Variant is the perturbedparameters that are a matrix (130x9)
>
> n = length(Paramvect);
> Variant = zeros((n^3+5*n+6)/6,n);
> p = n:1:1;
> s = 0;
> for k = n:1:n3 % Up to triples
> c = nchoosek(p,k);
> for m = 1:size(c,1)
> s = s + 1;
> Variant(s,c(m,:)) = Paramvect (c(m,:));
> end
> end
>
> I need to use this matrix that you called “Variant” and then run sensitivity analysis for each row of this matrix to get 130 different datasets.
>
> % Sensitivity analysis
> PcarConf.SolverOptions.SensitivityAnalysis = true;
> set (PcarConf.SensitivityAnalysisOptions, 'SpeciesOutputs', ...
> sbioselect(Pcar, 'Type', 'species'));
> get (PcarConf.SensitivityAnalysisOptions, 'SpeciesOutputs')
> set(PcarConf.SensitivityAnalysisOptions,'SpeciesInputFactors', ...
> sbioselect(Pcar, 'Type', 'species'));
> get (PcarConf.SensitivityAnalysisOptions, 'SpeciesInputFactors')
> set(PcarConf.SensitivityAnalysisOptions,'ParameterInputFactors', ...
> sbioselect(Pcar, 'Type', 'Parameter'));
> get (PcarConf.SensitivityAnalysisOptions, 'ParameterInputFactors')
>
> sensopts = PcarConf.SensitivityAnalysisOptions;
> paramvect = sbioselect(Pcar, 'Type', 'parameter');
> sensopts.ParameterInputFactors = paramvect;
> sensopts.Normalization = 'Full';
>
> for g=1:size(Variant,1)
> sdSens (g) = sbiosimulate(Pcar(g));
> [t(g),r(g),out(g),in(g)] = getsensmatrix(sdSens);
> end
>
> Thanks,
> Sam
