master

分支 (1)

管理

管理

master

sa-dl-code
/
MORRIS_GLOBALb_CFR_DL.m

close all;
clear;
clc;

rng('default');

krek2_origin=5000/10;
kBA1_origin=175;
kBA2_origin=35;
kbB_origin=250;
kbLF_origin=10;
kHCL_origin=7;
kbHC_origin=10;
kbF_origin=10;
kbct_origin=100;
kbNDH_origin=1;
krek1_origin=10000/10;
kbR_origin=10;
kbX_origin=10;
kbFd_origin=5;

percent=0.2;

krek2_range=[krek2_origin*(1-percent) krek2_origin*(1+percent)];
kBA1_range=[kBA1_origin*(1-percent) kBA1_origin*(1+percent)];
kBA2_range=[kBA2_origin*(1-percent) kBA2_origin*(1+percent)];
kbB_range=[kbB_origin*(1-percent) kbB_origin*(1+percent)];
kbLF_range=[kbLF_origin*(1-percent) kbLF_origin*(1+percent)];
kHCL_range=[kHCL_origin*(1-percent) kHCL_origin*(1+percent)];
kbHC_range=[kbHC_origin*(1-percent) kbHC_origin*(1+percent)];
kbF_range=[kbF_origin*(1-percent) kbF_origin*(1+percent)];
kbct_range=[kbct_origin*(1-percent) kbct_origin*(1+percent)];
kbNDH_range=[kbNDH_origin*(1-percent) kbNDH_origin*(1+percent)];
krek1_range=[krek1_origin*(1-percent) krek1_origin*(1+percent)];
kbR_range=[kbR_origin*(1-percent) kbR_origin*(1+percent)];
kbX_range=[kbX_origin*(1-percent) kbX_origin*(1+percent)];
kbFd_range=[kbFd_origin*(1-percent) kbFd_origin*(1+percent)];

PAB=1;
PAmB=0;
PABm=0;
PAmBm=0;
PABmm=0;
PAmBmm=0;
PpAB=0;
PpAmB=0;
PpABm=0;
PpAmBm=0;
PpABmm=0;
PpAmBmm=0;
S0=0.25;
S1=0.75;
S2=0;
S3=0;
PQ=2.5;
PQH=2.5;
LHF=1;
LmHF=0;
LHmF=0;
LHFm=0;
LmHmF=0;
LmHFm=0;
LHmFm=0;
LHmmF=0;
LmHmFm=0;
LmHmmF=0;
LHmmFm=0;
LmHmmFm=0;
PC=3;
PCp=0;
Fd=3;
Fdm=0;
RX=0.63;
RXm=0;
RpX=0;
RpXm=0;
FNRi=3;
FNRa=0;
FNRam=0;
FNRamm=0;
PQm=0;

x0=[PAB;PAmB;PABm;PAmBm;PABmm;PAmBmm;...
    PpAB;PpAmB;PpABm;PpAmBm;PpABmm;PpAmBmm;...
    S0;S1;S2;S3;...
    PQ;PQH;...
    LHF;LmHF;LHmF;LHFm;LmHmF;LmHFm;LHmFm;LHmmF;LmHmFm;LmHmmF;LHmmFm;LmHmmFm;...
    PC;PCp;...
    Fd;Fdm;...
    RX;RXm;RpX;RpXm;...
    FNRi;FNRa;FNRam;FNRamm;...
    PQm];

%% Perform Morris sensitivity analysis

num_levels = 4;
delta = num_levels / (2*(num_levels -1));

krek2_samples = linspace(0, 1, num_levels)';
kBA1_samples = linspace(0, 1, num_levels)';
kBA2_samples = linspace(0, 1, num_levels)';
kbB_samples = linspace(0, 1, num_levels)';
kbLF_samples = linspace(0, 1, num_levels)';
kHCL_samples = linspace(0, 1, num_levels)';
kbHC_samples = linspace(0, 1, num_levels)';
kbF_samples = linspace(0, 1, num_levels)';
kbct_samples = linspace(0, 1, num_levels)';
kbNDH_samples = linspace(0, 1, num_levels)';
krek1_samples = linspace(0, 1, num_levels)';
kbR_samples = linspace(0, 1, num_levels)';
kbX_samples = linspace(0, 1, num_levels)';
kbFd_samples = linspace(0, 1, num_levels)';

num_parameters = 14;

parameter_samples_matrix = [krek2_samples, kBA1_samples, kBA2_samples, kbB_samples, kbLF_samples, kHCL_samples, kbHC_samples, kbF_samples, kbct_samples, kbNDH_samples, krek1_samples, kbR_samples, kbX_samples, kbFd_samples];

x_star_matrix = zeros(num_levels, num_parameters); % X* initial matrix

for i = 1:num_parameters
    x_star_matrix(:, i) = datasample(parameter_samples_matrix(:, i), num_levels, 'Replace', false);
end

B_matrix = tril(ones(15, 14), -1);

J_matrix_15times14 = ones(15, 14);
J_matrix_15times1 = ones(15, 1);

D_star_matrix = diag(2 * randi([0, 1], 14, 1) - 1);
% sequence = ones(1, 18);
% sequence(2:2:end) = -1;
% D_star_matrix = diag(sequence);

P_star_matrix = eye(14);
% P_star_matrix = P_star(:, randperm(18));

B_star_matrix = (J_matrix_15times1*x_star_matrix(4,:)+delta/2*((2*B_matrix-J_matrix_15times14)*D_star_matrix+J_matrix_15times14))*P_star_matrix;

num_Bstar_parameters = size(B_star_matrix, 2); %  14
num_Bstar_samples = size(B_star_matrix, 1); %  15

mapped_values_matrix = zeros(num_Bstar_samples, num_Bstar_parameters);

param_ranges = {krek2_range, kBA1_range, kBA2_range, kbB_range, kbLF_range, kHCL_range, kbHC_range, kbF_range, kbct_range, kbNDH_range, krek1_range, kbR_range, kbX_range, kbFd_range};

for i = 1:num_Bstar_parameters
    mapped_values_matrix(:, i) = param_ranges{i}(1) + B_star_matrix(:, i) * (param_ranges{i}(2) - param_ranges{i}(1));
end

t=0:1e-4:1;

qyFtotal_values = zeros(num_Bstar_samples, length(t));

for i = 1:num_Bstar_samples

        [t, x] = ode15s(@(t, x) ODE_MORRIS_GLOBALb_CFR_DL(t, x, mapped_values_matrix(i,:)), t, x0);

        PSIIm = x(:,2) + x(:,4) + x(:,6) + x(:,8) + x(:,10) + x(:,12);
        PSIIo = x(:,1) + x(:,3) + x(:,5) + x(:,7) + x(:,9) + x(:,11);
        PSIm = x(:,36) + x(:,38);
        qyF_PSIIopen = 0.02;
        qyF_PSIIclosed = 0.08;
        qyFPSII = qyF_PSIIopen * PSIIo + qyF_PSIIclosed * PSIIm;
        qyF_PSIclosed = 0.15 * qyF_PSIIclosed;
        qyFPSI = qyF_PSIclosed * PSIm;
        qyFtotal = qyFPSII + qyFPSI;

        qyFtotal_values(i, :) = qyFtotal;
end

percentiles = prctile(qyFtotal_values, [5, 50, 95]);

%% Calculate Morris Indices

EE_krek2_matrix = zeros(length(t), 1);
EE_kBA1_matrix = zeros(length(t), 1);
EE_kBA2_matrix = zeros(length(t), 1);
EE_kbB_matrix = zeros(length(t), 1);
EE_kbLF_matrix = zeros(length(t), 1);
EE_kHCL_matrix = zeros(length(t), 1);
EE_kbHC_matrix = zeros(length(t), 1);
EE_kbF_matrix = zeros(length(t), 1);
EE_kbct_matrix = zeros(length(t), 1);
EE_kbNDH_matrix = zeros(length(t), 1);
EE_krek1_matrix = zeros(length(t), 1);
EE_kbR_matrix = zeros(length(t), 1);
EE_kbX_matrix = zeros(length(t), 1);
EE_kbFd_matrix = zeros(length(t), 1);

if B_star_matrix(2, 1) > B_star_matrix(1, 1)
    EE_krek2_matrix(:, 1) = (qyFtotal_values(2, :) - qyFtotal_values(1, :)) / delta;
else
    EE_krek2_matrix(:, 1) = (qyFtotal_values(1, :) - qyFtotal_values(2, :)) / delta;
end

if B_star_matrix(3, 2) > B_star_matrix(2, 2)
    EE_kBA1_matrix(:, 1) = (qyFtotal_values(3, :) - qyFtotal_values(2, :)) / delta;
else
    EE_kBA1_matrix(:, 1) = (qyFtotal_values(2, :) - qyFtotal_values(3, :)) / delta;
end

if B_star_matrix(4, 3) > B_star_matrix(3, 3)
    EE_kBA2_matrix(:, 1) = (qyFtotal_values(4, :) - qyFtotal_values(3, :)) / delta;
else
    EE_kBA2_matrix(:, 1) = (qyFtotal_values(3, :) - qyFtotal_values(4, :)) / delta;
end

if B_star_matrix(5, 4) > B_star_matrix(4, 4)
    EE_kbB_matrix(:, 1) = (qyFtotal_values(5, :) - qyFtotal_values(4, :)) / delta;
else
    EE_kbB_matrix(:, 1) = (qyFtotal_values(4, :) - qyFtotal_values(5, :)) / delta;
end

if B_star_matrix(6, 5) > B_star_matrix(5, 5)
    EE_kbLF_matrix(:, 1) = (qyFtotal_values(6, :) - qyFtotal_values(5, :)) / delta;
else
    EE_kbLF_matrix(:, 1) = (qyFtotal_values(5, :) - qyFtotal_values(6, :)) / delta;
end

if B_star_matrix(7, 6) > B_star_matrix(6, 6)
    EE_kHCL_matrix(:, 1) = (qyFtotal_values(7, :) - qyFtotal_values(6, :)) / delta;
else
    EE_kHCL_matrix(:, 1) = (qyFtotal_values(6, :) - qyFtotal_values(7, :)) / delta;
end

if B_star_matrix(8, 7) > B_star_matrix(7, 7)
    EE_kbHC_matrix(:, 1) = (qyFtotal_values(8, :) - qyFtotal_values(7, :)) / delta;
else
    EE_kbHC_matrix(:, 1) = (qyFtotal_values(7, :) - qyFtotal_values(8, :)) / delta;
end

if B_star_matrix(9, 8) > B_star_matrix(8, 8)
    EE_kbF_matrix(:, 1) = (qyFtotal_values(9, :) - qyFtotal_values(8, :)) / delta;
else
    EE_kbF_matrix(:, 1) = (qyFtotal_values(8, :) - qyFtotal_values(9, :)) / delta;
end

if B_star_matrix(10, 9) > B_star_matrix(9, 9)
    EE_kbct_matrix(:, 1) = (qyFtotal_values(10, :) - qyFtotal_values(9, :)) / delta;
else
    EE_kbct_matrix(:, 1) = (qyFtotal_values(9, :) - qyFtotal_values(10, :)) / delta;
end

if B_star_matrix(11, 10) > B_star_matrix(10, 10)
    EE_kbNDH_matrix(:, 1) = (qyFtotal_values(11, :) - qyFtotal_values(10, :)) / delta;
else
    EE_kbNDH_matrix(:, 1) = (qyFtotal_values(10, :) - qyFtotal_values(11, :)) / delta;
end

if B_star_matrix(12, 11) > B_star_matrix(11, 11)
    EE_krek1_matrix(:, 1) = (qyFtotal_values(12, :) - qyFtotal_values(11, :)) / delta;
else
    EE_krek1_matrix(:, 1) = (qyFtotal_values(11, :) - qyFtotal_values(12, :)) / delta;
end

if B_star_matrix(13, 12) > B_star_matrix(12, 12)
    EE_kbR_matrix(:, 1) = (qyFtotal_values(13, :) - qyFtotal_values(12, :)) / delta;
else
    EE_kbR_matrix(:, 1) = (qyFtotal_values(12, :) - qyFtotal_values(13, :)) / delta;
end

if B_star_matrix(14, 13) > B_star_matrix(13, 13)
    EE_kbX_matrix(:, 1) = (qyFtotal_values(14, :) - qyFtotal_values(13, :)) / delta;
else
    EE_kbX_matrix(:, 1) = (qyFtotal_values(13, :) - qyFtotal_values(14, :)) / delta;
end

if B_star_matrix(15, 14) > B_star_matrix(14, 14)
    EE_kbFd_matrix(:, 1) = (qyFtotal_values(15, :) - qyFtotal_values(14, :)) / delta;
else
    EE_kbFd_matrix(:, 1) = (qyFtotal_values(14, :) - qyFtotal_values(15, :)) / delta;
end

%% Plot the results

figure(1);
hold on;

for i = 1:num_Bstar_samples
    plot(t, qyFtotal_values(i, :), 'Color', [0.7, 0.7, 0.7]);
end

median_line = plot(t, percentiles(2, :), 'b', 'LineWidth', 1.5);

quartile_lines_90 = plot(t, [percentiles(1, :); percentiles(3, :)], 'r', 'LineWidth', 1.5);

set(gca, 'XScale', 'log');
xlabel('Time (s)');
ylabel('Simulated Fluorescence Emission (a.u.)');
hold off;

legend([median_line, quartile_lines_90(1)], {'Median (50%)', '5%-95% Range'}, 'Location', 'NorthWest');

legend('boxoff');
hLegend = findobj(gcf, 'Type', 'Legend');
set(hLegend, 'FontSize', 10);

% resolution = 600;
% print('output_image', '-dpng', ['-r', num2str(resolution)]);