В этом примере показано, как вычислить чувствительность некоторых разновидностей в модели Лотки-Вольтерры с помощью SimFunctionSensitivity объект.

Загрузите демонстрационный проект.

sbioloadproject lotka;

Задайте входные параметры.

params = {'Reaction1.c1', 'Reaction2.c2'};

Задайте наблюдаемые разновидности, которые являются выходными параметрами симуляции.

observables  = {'y1', 'y2'};

Создайте SimFunctionSensitivity объект. Установите чувствительность выходные факторы на все разновидности (y1 и y2) заданный в observables аргумент и входные факторы тем в params аргумент (c1 и c2) путем установки аргумента пары "имя-значение" 'all'.

f = createSimFunction(m1,params,observables,[],'SensitivityOutputs','all','SensitivityInputs','all','SensitivityNormalization','Full')

f = 
SimFunction

Parameters:

          Name          Value        Type     
    ________________    _____    _____________

    {'Reaction1.c1'}      10     {'parameter'}
    {'Reaction2.c2'}    0.01     {'parameter'}

Observables: 

     Name        Type    
    ______    ___________

    {'y1'}    {'species'}
    {'y2'}    {'species'}

Dosed: None

Sensitivity Input Factors: 

          Name              Type     
    ________________    _____________

    {'Reaction1.c1'}    {'parameter'}
    {'Reaction2.c2'}    {'parameter'}

Sensitivity Output Factors: 

     Name        Type    
    ______    ___________

    {'y1'}    {'species'}
    {'y2'}    {'species'}

Sensitivity Normalization: 

Full

Вычислите чувствительность путем выполнения объекта с c1 и c2 установите на 10 и 0.1, соответственно. Установите выходные времена от 1 до 10. t содержит моменты времени, y содержит данные моделирования и sensMatrix матрица чувствительности, содержащая чувствительность y1 и y2 относительно c1 и c2.

[t,y,sensMatrix] = f([10,0.1],[],[],1:10);

Получите информацию о чувствительности в моменте времени 5.

temp = sensMatrix{:};
sensMatrix2 = temp(t{:}==5,:,:);
sensMatrix2 = squeeze(sensMatrix2)

sensMatrix2 = 2×2

   37.6987   -6.8447
  -40.2791    5.8225

Строки sensMatrix2 представляйте выходные факторы (y1 и y2). Столбцы представляют входные факторы (c1 и c2).

Установите время остановки на 15, не задавая выходные времена. В этом случае выходные времена являются моментами времени решателя по умолчанию.

sd = f([10,0.1],15);

Получите расчетную чувствительность из SimData объект sd.

[t,y,outputs,inputs] = getsensmatrix(sd);

Постройте чувствительность разновидностей y1 и y2 относительно c1.

figure;
plot(t,y(:,:,1));
legend(outputs);
title('Sensitivities of species y1 and y2 with respect to parameter c1');
xlabel('Time');
ylabel('Sensitivity');

Постройте чувствительность разновидностей y1 и y2 относительно c2.

figure;
plot(t,y(:,:,2));
legend(outputs);
title('Sensitivities of species y1 and y2 with respect to parameter c2');
xlabel('Time');
ylabel('Sensitivity');

В качестве альтернативы можно использовать sbioplot.

sbioplot(sd);

Можно также построить матрицу чувствительности использование интеграла времени для расчетной чувствительности y1 и y2. График показывает y1 и y2 более чувствительны к c1 чем c2.

[~, in, out] = size(y);
result = zeros(in, out);
for i = 1:in
    for j = 1:out
        result(i,j) = trapz(t(:),abs(y(:,i,j)));
    end
end
figure;
hbar = bar(result);
haxes = hbar(1).Parent;
haxes.XTick = 1:length(outputs);
haxes.XTickLabel = outputs;
legend(inputs,'Location','NorthEastOutside');
ylabel('Sensitivity');