What the phase C values are not coming Rpeak = 90% and Rss = 58% while other phase values are correct. Kindly explain is their any logical error in the code

What the phase C values are not coming Rpeak = 90% and Rss = 58% while other phase values are correct. Kindly explain is their any logical error in the code. Thanks
% Define parameters
Kf_Max = 100; % maximum forward rate
RLC = [0.01, 1, 10,0.01]; % RLC values
TauKf_ON = -0.01; % TauKf_ON
TauKf_OFF = -0.01; % TauKf_OFF
Kb_Max = 80; % maximum backward rate
Kb_Min = 10; % minimum backward rate
TauKb_ON = -0.01; % TauKb_ON
TauKb_OFF = -0.01; % TauKb_OFF
PhaseTimes = [0, 10, 500, 1000, 2000];% phase times (each row defines a phase)
timespan = [0, 2000]; % timespan for the simulation

% Example initial conditions for non-active and active receptor concentrations
Non_Active_Receptor_concentration = 100; % Initial non-active receptor concentration (as a vector)
Active_Receptor_concentration =0; % Initial active receptor concentration (as a vector)

[t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] =TestCode1(Kf_Max, RLC, TauKf_ON, …
Kb_Max, Kb_Min, TauKb_ON, PhaseTimes, timespan, Non_Active_Receptor_concentration, Active_Receptor_concentration);

function [t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] =TestCode1(Kf_Max, RLC, TauKf_ON, …
Kb_Max, Kb_Min, TauKb_ON, PhaseTimes, timespan, Non_Active_Receptor_concentration, Active_Receptor_concentration)

% Define functions for forward and backward reaction rates
Kf_L = @(t) calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON);
Kb = @(t) calculate_kb(t, PhaseTimes, Kb_Max, Kb_Min, TauKb_ON, RLC);

% Ensure that Non_Active_Receptor_concentration and Active_Receptor_concentration are column vectors
Non_Active_Receptor_concentration = Non_Active_Receptor_concentration(:);
Active_Receptor_concentration = Active_Receptor_concentration(:);

% Set initial conditions
initial_conditions = [Non_Active_Receptor_concentration(1); Active_Receptor_concentration(1)];

% Set ODE options for step size
options = odeset(‘MaxStep’, 0.05, ‘RelTol’, 1e-6, ‘AbsTol’, 1e-8);

% Solve the ODE system
[t, y] = ode45(@(t, y) ode_LR(t, y, Kf_L, Kb), timespan, initial_conditions, options);

% Extract the concentrations
Non_Active_Receptor_concentration = y(:, 1);
Active_Receptor_concentration = y(:, 2);

% Calculate forward and backward reaction rates over time
Kf_values = arrayfun(Kf_L, t);
Kb_values = arrayfun(Kb, t);

% Plot active and non-active receptor concentrations
figure;
plot(t, Non_Active_Receptor_concentration, ‘r’, ‘DisplayName’, ‘Non-Active Receptor Concentration’);
hold on;
plot(t, Active_Receptor_concentration, ‘b’, ‘DisplayName’, ‘Active Receptor Concentration’);
legend;
xlabel(‘Time’);
ylabel(‘Concentration’);
title(‘Receptor Concentrations’);
hold off;

% Plot forward and backward reaction rates
figure;
plot(t, Kf_values, ‘k’, ‘DisplayName’, ‘Forward Reaction Rate’);
hold on;
plot(t, Kb_values, ‘c’, ‘DisplayName’, ‘Backward Reaction Rate’);
legend;
xlabel(‘Time’);
ylabel(‘Reaction Rate’);
title(‘Reaction Rates’);
hold off;

% Calculate phase results
PhaseResults = arrayfun(@(i) …
calculate_phase(t, Active_Receptor_concentration, PhaseTimes(:,i:i+1), RLC(i)) …
, 1:(size(PhaseTimes, 2)-1), ‘UniformOutput’, false);
PhaseResults = vertcat(PhaseResults{:}); % Concatenate results

% Calculate maximum forward reaction rate
Kf_LMax = Kf_Max * (RLC ./ (RLC + 1));

% Write Phase results to CSV
for i = 1:size(PhaseResults, 1)
PhaseTable = struct2table(PhaseResults(i))
writetable(PhaseTable, [‘Phase’, num2str(i), ‘.csv’]);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kf_L = calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON)
Kf_LMax = Kf_Max * (RLC ./ (RLC + 1));
Kf_L = Kf_LMax(1); % Default to the first phase value

num_phases = numel(RLC);
for j = 1:num_phases
T_start = PhaseTimes(j);
T_end = PhaseTimes(j + 1);
if t >= T_start && t < T_end
if j == 1
Kf_L = Kf_LMax(j);
else
prev_end = PhaseTimes(j);
if j < num_phases
if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
Kf_L = Kf_LMax(j);
elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
Kf_end = Kf_LMax(j) – (Kf_LMax(j) – Kf_LMax(j-1)) * exp(TauKf_ON * (T_end – T_start));
Kf_L = Kf_LMax(j) + (Kf_end – Kf_LMax(j-1)) * exp(TauKf_ON * (t – prev_end));
elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
Kf_end = Kf_LMax(j) – (Kf_LMax(j) – Kf_LMax(j-1)) * exp(TauKf_ON * (T_end – T_start));
Kf_L = Kf_LMax(j) + (Kf_end – Kf_LMax(j-1)) * exp(TauKf_ON * (t – prev_end));
elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
Kf_end = Kf_LMax(j) – (Kf_LMax(j) – Kf_LMax(j-1)) * exp(TauKf_ON * (T_end – T_start));
Kf_L = Kf_LMax(j) + (Kf_end – Kf_LMax(j-1)) * exp(TauKf_ON * (t – prev_end));
end
end
end
return;
end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kb = calculate_kb(t, PhaseTimes, Kb_Max, Kb_Min, ~, RLC)
Kb = Kb_Max; % Default to the minimum value
if all(RLC == RLC(1))
Kb = Kb_Min; % Set Kb to Kb_Min if all RLC values are equal
return;
end
for j = 1:numel(RLC)
T_start = PhaseTimes(j);
T_end = PhaseTimes(j + 1);
if t >= T_start && t < T_end
if j == 1
Kb = Kb_Max;
else
% prev_end = PhaseTimes(j);
if j < numel(RLC)
if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
Kb = Kb_Max ;
elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
Kb = Kb_Max ;
elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
Kb = Kb_Max;
elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
Kb = Kb_Max ;
end
else
if RLC(j-1) < RLC(j)
Kb = Kb_Max;

elseif RLC(j-1) > RLC(j)
Kb = Kb_Max;
end
end
end
return;
end
end
end

dNonActiveReceptor_dt = -Kf_L(t) * Non_Active_Receptor_concentration + Kb(t) * Active_Receptor_concentration;
dActiveReceptor_dt = Kf_L(t) * Non_Active_Receptor_concentration – Kb(t) * Active_Receptor_concentration;

dydt = [dNonActiveReceptor_dt; dActiveReceptor_dt];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result = calculate_phase(t, Active_Receptor_concentration, PhaseTimes, RLC)
T_start = PhaseTimes(:,1);
T_end = PhaseTimes(:,2);
Phase_indices = (t >= T_start & t <= T_end);
Phase_concentration = Active_Receptor_concentration(Phase_indices);
Phase_time = t(Phase_indices);

% Calculate peak and steady-state values
[Rpeak, Tpeak, peak_type] = findpeak(Phase_time, Phase_concentration);
Rss = mean(Active_Receptor_concentration(t >= (T_end – 10) & t <= T_end));

% Calculate the T50 (time to reach half of peak value)
half_peak_value = (Rss + Rpeak) / 2;
[~, idx_50_percent] = min(abs(Phase_concentration – half_peak_value));
T50 = Phase_time(idx_50_percent) – Tpeak;

% Calculate other metrics
ratio_Rpeak_Rss = Rpeak / Rss;
Delta = Rpeak – Rss;
L = RLC;

% Package results in a struct
result.Rpeak = Rpeak;
result.Rss = Rss;
result.Tpeak = Tpeak;
result.T50 = T50;
result.ratio_Rpeak_Rss = ratio_Rpeak_Rss;
result.Delta = Delta;
result.L = L;
result.peak_type = peak_type;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Rpeak, Tpeak, peak_type] = findpeak(time, concentration)
% Compute the derivative
dt = diff(time);
dy = diff(concentration);
derivative = dy ./ dt;

% Find zero-crossings of the derivative
zero_crossings = find(diff(sign(derivative)) ~= 0);

% Initialize output variables
Rpeak = NaN;
Tpeak = NaN;
peak_type = ‘None’;

% Check if there are zero crossings
if ~isempty(zero_crossings)
% Identify peaks and troughs by examining the sign changes
max_indices = zero_crossings(derivative(zero_crossings) > 0 & derivative(zero_crossings + 1) < 0);
min_indices = zero_crossings(derivative(zero_crossings) < 0 & derivative(zero_crossings + 1) > 0);

% Check if there are maxima or minima
if ~isempty(max_indices) || ~isempty(min_indices)
if ~isempty(max_indices) && ~isempty(min_indices)
% Find the highest maximum
[Rpeak_max, maxIdx] = max(concentration(max_indices));
% Find the lowest minimum
[Rpeak_min, minIdx] = min(concentration(min_indices));

% Determine whether the highest maximum is greater or the lowest minimum is smaller
if Rpeak_max >= abs(Rpeak_min)
Rpeak = Rpeak_max;
Tpeak = time(max_indices(maxIdx));
peak_type = ‘High’;
else
Rpeak = Rpeak_min;
Tpeak = time(min_indices(minIdx));
peak_type = ‘Low’;
end
% If only maxima exist
elseif ~isempty(max_indices)
[Rpeak, maxIdx] = max(concentration(max_indices));
Tpeak = time(max_indices(maxIdx));
peak_type = ‘High’;
% If only minima exist
elseif ~isempty(min_indices)
[Rpeak, minIdx] = min(concentration(min_indices));
Tpeak = time(min_indices(minIdx));
peak_type = ‘Low’;
end
end
end
endWhat the phase C values are not coming Rpeak = 90% and Rss = 58% while other phase values are correct. Kindly explain is their any logical error in the code. Thanks
% Define parameters
Kf_Max = 100; % maximum forward rate
RLC = [0.01, 1, 10,0.01]; % RLC values
TauKf_ON = -0.01; % TauKf_ON
TauKf_OFF = -0.01; % TauKf_OFF
Kb_Max = 80; % maximum backward rate
Kb_Min = 10; % minimum backward rate
TauKb_ON = -0.01; % TauKb_ON
TauKb_OFF = -0.01; % TauKb_OFF
PhaseTimes = [0, 10, 500, 1000, 2000];% phase times (each row defines a phase)
timespan = [0, 2000]; % timespan for the simulation