function [t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = TestCode1(Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...
Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, 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, TauKf_OFF);
Kb = @(t) calculate_kb(t, PhaseTimes, Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, RLC); % Pass RLC here
% 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;
% Extract values for each phase
if size(PhaseTimes, 1) ~= numel(RLC)
error('PhaseTimes and RLC must have the same number of elements.');
end
PhaseResults = arrayfun(@(i) calculate_phase(t, Active_Receptor_concentration, PhaseTimes(i, :), RLC(i)), 1:size(PhaseTimes, 1), 'UniformOutput', false);
PhaseResults = vertcat(PhaseResults{:});
Kf_LMax = Kf_Max * (RLC ./ (RLC + 1));
% Plot peak values
figure;
plot(t, Active_Receptor_concentration, 'm', 'DisplayName', 'Active Receptor Concentration');
hold on;
% Convert PhaseResults to a struct array if it is a cell array
if iscell(PhaseResults)
PhaseResults = cell2mat(PhaseResults);
end
for i = 1:size(PhaseResults, 1)
plot(PhaseResults(i).Tpeak, PhaseResults(i).Rpeak, 'o', 'DisplayName', ['Peak ', num2str(i)]);
end
legend;
xlabel('Time');
ylabel('Concentration');
title('Active Receptor Concentration with Peaks');
hold off;
% 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 = size(PhaseTimes, 1);
for j = 1:num_phases
T_start = PhaseTimes(j, 1);
T_end = PhaseTimes(j, 2);
if t >= T_start && t < T_end
if j == 1
Kf_L = Kf_LMax(j);
else
prev_end = PhaseTimes(j-1, 2);
if j < length(RLC)
% Ensure indices j+1 are within bounds
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_L = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
% Kf_L = Kf_LMax(j) - (Kf_endA - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
% elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
% Kf_L = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
% elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
% Kf_L = Kf_LMax(j) - (Kf_LMax(j) - 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, TauKb_ON, ~, RLC)
Kb = Kb_Min; % Default to the minimum value
% Check if all RLC values are equal
if all(RLC == RLC(1))
Kb = Kb_Min; % Set Kb to Kb_Min if all RLC values are equal
return;
end
for j = 1:size(PhaseTimes, 1)
T_start = PhaseTimes(j, 1);
T_end = PhaseTimes(j, 2);
if t >= T_start && t < T_end
if j == 1
Kb = Kb_Min;
else
prev_end = PhaseTimes(j-1, 2);
% Add conditions for smooth behavior based on RLC patterns
if j < length(RLC)
if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
% RLC increases then decreases
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
% RLC increases then continues to increase
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
% RLC decreases then increases
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end)); % Change start from Kb_Max
elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
% RLC decreases then continues to decrease
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end)); % Change start from Kb_Max
end
else
if RLC(j-1) < RLC(j)
% RLC increases
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
elseif RLC(j-1) > RLC(j)
% RLC decreases
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end)); % Change start from Kb_Max
end
end
end
return;
end
end
end
% function Kb = calculate_kb(t, PhaseTimes, Kb_Max, Kb_Min, TauKb_ON, ~, RLC)
% Kb = Kb_Min; % Default to the minimum value
% % Check if all RLC values are equal
% if all(RLC == RLC(1))
% Kb = Kb_Min; % Set Kb to Kb_Min if all RLC values are equal
% return;
% end
% for j = 1:size(PhaseTimes, 1)
% T_start = PhaseTimes(j, 1);
% T_end = PhaseTimes(j, 2);
% if t >= T_start && t < T_end
% if j == 1
% Kb = Kb_Min;
% else
% prev_end = PhaseTimes(j-1, 2);
% % Add conditions for smooth behavior based on RLC patterns
% if j < length(RLC)
% if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
% % RLC increases then decreases
% Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
% elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
% % RLC increases then continues to increase
% Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
% elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
% % RLC decreases then increases
% Kb = Kb_Min + (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
% elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
% % RLC decreases then continues to decrease
% Kb = Kb_Min + (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
% end
% else
% if RLC(j-1) < RLC(j)
% % RLC increases
% Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
% elseif RLC(j-1) > RLC(j)
% % RLC decreases
% Kb = Kb_Min + (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
% end
% end
% end
% return;
% end
% 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 = size(PhaseTimes, 1);
% for j = 1:num_phases
% T_start = PhaseTimes(j, 1);
% T_end = PhaseTimes(j, 2);
% if t >= T_start && t < T_end
% if j == 1
% Kf_L = Kf_LMax(j);
% else
% prev_end = PhaseTimes(j-1, 2);
% if j < num_phases
% % Ensure indices j+1 are within bounds
% if j + 1 <= length(RLC) && RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
% Kf_L = Kf_LMax(j);
% elseif j + 1 <= length(RLC) && RLC(j-1) < RLC(j) && RLC(j) < RLC(j+1)
% Kf_L = Kf_LMax(j) - (Kf_LMax(j) - 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, TauKb_ON, ~, RLC)
% Kb = Kb_Min; % Default to the minimum value
% % Check if all RLC values are equal
% if all(RLC == RLC(1))
% Kb = Kb_Min; % Set Kb to Kb_Min if all RLC values are equal
% return;
% end
% for j = 1:size(PhaseTimes, 1)
% T_start = PhaseTimes(j, 1);
% T_end = PhaseTimes(j, 2);
% if t >= T_start && t < T_end
% if j == 1
% Kb = Kb_Min;
% else
% prev_end = PhaseTimes(j-1, 2);
% % Add conditions for smooth behavior based on RLC patterns
% if j + 1 <= length(RLC) && RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
% Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
% elseif j + 1 <= length(RLC) && RLC(j-1) < RLC(j) && RLC(j) < RLC(j+1)
% Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
% end
% end
% return;
% end
% end
% end
function dydt = ode_LR(t, y, Kf_L, Kb)
Non_Active_Receptor_concentration = y(1);
Active_Receptor_concentration = y(2);
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
end
