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

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

sbioloadproject lotka;

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

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

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

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

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

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. Расширьте Run1, чтобы выбрать который симуляция или данные о чувствительности, чтобы отобразиться.

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