Sensitivity analysis problem in simple prey-predator model

34 views Asked by At

I am trying to perform sensitivity analysis by using PRCC and latin hypercube sampling. I wanted to take a general model, so I chose the prey-predator. Mathematically, I used ODEs system which is simulated by ode45. However, I am getting an empty plot for PRCC valus against the sampled parameters (I want to test them all). What is wrong in this approach? how do I change it?



% Prey-predator ODEs model with sensitivity analysis and plotting of statistically significant PRCC values

    % Time span
    tspan = [0 10];

    % Initial conditions
    y0 = [40; 9];   % Prey initial count, Predator initial count

    % Latin hypercube sampling
    numSamples = 1000;
    parameters = lhsdesign(numSamples, 4);

    % Define parameter ranges
    alphaRange = linspace(0.05, 0.15, numSamples)';
    betaRange = linspace(0.01, 0.03, numSamples)';
    deltaRange = linspace(0.01, 0.03, numSamples)';
    gammaRange = linspace(0.1, 0.3, numSamples)';

    % Scale the Latin hypercube samples to the parameter ranges
    alpha = interp1(linspace(0, 1, numSamples), alphaRange, parameters(:, 1));
    beta = interp1(linspace(0, 1, numSamples), betaRange, parameters(:, 2));
    delta = interp1(linspace(0, 1, numSamples), deltaRange, parameters(:, 3));
    gamma = interp1(linspace(0, 1, numSamples), gammaRange, parameters(:, 4));
    
    parameters = [alpha, beta, delta, gamma];

    % Sensitivity analysis
    preyCountAtDay10 = zeros(numSamples, 1);
    for i = 1:numSamples
        [t, y] = ode45(@(t, y) odeSystem(t, y, parameters(i, :)), tspan, y0);
        preyCountAtDay10(i) = y(end, 1); % Prey count at day 10
    end

    % PRCC calculation
    [prcc, pValues] = partialRankCorrelationCoefficient(preyCountAtDay10, parameters(:, 1));

    % Find statistically significant PRCC values (p-value < 0.01)
    significantPRCC = prcc(pValues < 0.01);
    significantParameters = parameters(pValues < 0.01, :);
    
    % Create a bar plot of statistically significant PRCC values
    figure;
    bar(significantPRCC);
    ylabel('PRCC for Prey Count at Day 10');
    xlabel('Indices of Statistically Significant Parameters');
    title('Statistically Significant PRCC Values');
    xticklabels({'Parameter 1', 'Parameter 2', 'Parameter 3', 'Parameter 4'});

% ODE system
function dydt = odeSystem(t, y, parameters)
    alpha = parameters(1);
    beta = parameters(2);
    delta = parameters(3);
    gamma = parameters(4);

    prey = y(1);
    predator = y(2);

    preyGrowth = alpha * prey - beta * prey * predator;
    predatorGrowth = gamma * prey * predator - delta * predator;

    dydt = [preyGrowth; predatorGrowth];
end

% Partial Rank Correlation Coefficient Calculation
function [prcc, pValues] = partialRankCorrelationCoefficient(y, x)
    % Calculate Pearson correlation coefficient
    r = corr(y, x);

    % Sample size
    n = length(y);

    % Calculate the t-statistic
    t = r * sqrt((n - 2) / (1 - r^2));

    % Calculate the p-value
    pValues = 2 * (1 - tcdf(abs(t), n - 2));

    % Calculate PRCC
    prcc = r * sqrt((n - 2) / (1 - r^2));
end

0

There are 0 answers