How i modified the code so Kb start form Kb_min for first phase but after that it start at Kb_Max . Kindly guide me

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