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