В этом примере показано, как рассчитать чувствительность некоторых видов в модели Lotka-Volterra с помощью 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');