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

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

type TestCode1.m

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

[Sahas]

Hi @Ehtisham,

Could you please provide more details about the requirements for this code or a reference graph related to it? This information would be very helpful for debugging purposes.

[Ehtisham]

i am getting this graph but i want that Kb in first phase just start form Kb_min but in other phase kept the Kb_Max like this graph blue line graph and also am getting the Just Phase1 Csv file not for all phase @Sahas

Answer 1

[Ketan]

you need to adjust the "calculate_kb" function accordingly