When i call and run this code it just save the Phase1 results in CSvV file not other results .

When i call and run this code it just save the Phase1 results in CSvV file not other results .

[Jaimin]

Hi @Ehtisham,

Could you please provide detailed information about the expected results in CSV format? Thank you.

[Ehtisham]

I need a csv file of each phase with following columns (Rpeak, Rss, Tpeak, T50 ratio_Rpeak_Rss, Delta, L, peak_type)

@Jaimin

Answer 1

[nick]

Hi Ehtisham,

I understand that you are unable to save CSV results other than 'Phase1.csv'.

The issue occurs because the PhaseResults variable is a 1x1 struct. This happened as the array used in 'arrayfun', as shown in the following code from TestCode1.m, contains only 1 element :

  % Calculate phase results
    PhaseResults = arrayfun(@(i) ...
        calculate_phase(t, Active_Receptor_concentration, PhaseTimes(i, :), RLC(i)) ...
        , 1:size(PhaseTimes, 1), 'UniformOutput', false);

The size of PhaseTimes, a row vector, is [1 8] which leads to the array utilized in the 'arrayfun' to have a single element. To resolve the issue, kindly ensure that the array of the correct dimension is created.

You may refer to the following documentation to learn more about 'arrayfun' function:

https://www.mathworks.com/help/matlab/ref/arrayfun.html

Hope this helps!

[Ehtisham]

% Calculate phase results

PhaseResults = arrayfun(@(i) calculate_phase(t, Active_Receptor_concentration, PhaseTimes(i, :), RLC(i)), 1:size(PhaseTimes), '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

now i write that command but again getting the same Phase1 csv file. @Neelanshu

[nick]

Hello @Ehtisham,

I understand that you are getting only Phase1.csv.

This occured as even after the changes the size of array used in the function 'arrayfun' contains 1 element that is 1. I have attached the image below of the output array for your reference.

To create an array with more than element with PhaseTimes, you may use the following code :

1:size(PhaseTimes, 2)

Hope this helps!

[Ehtisham]

I did not get you where i need to do that change. can you kindly elborate this @Neelanshu see it giving error

[nick]

Hi @Ehtisham,

Sure I will elaborate the cause of error. The issues occurs because the 'calculate_phase' function is not compatible with row vectors like 'PhaseTimes'.

I have updated the code at the following sections :

PhaseResults = arrayfun(@(i) ...

calculate_phase(t, Active_Receptor_concentration, PhaseTimes(:,i:i+1), RLC(i)) ...

, 1:(size(PhaseTimes, 2)-1), 'UniformOutput', false);

This sections creates a 1x7 array and calculate_phase works on the elements of this array.

function result = calculate_phase(t, Active_Receptor_concentration, PhaseTimes, RLC)

T_start = PhaseTimes(:,1);

T_end = PhaseTimes(:,2);

In this section, calculate_phase is modified to work with row vectors.

This results in the formation of 7 CSV files in the current path as shown in the image :

Here is the entire code script:

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.01, 1, 5, 1, 7, 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, 3000, 4000,5000];% phase times (each row defines a phase)
timespan = [0, 5000];  % 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)

% Run the 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);

PhaseTable = 1x8 table
    Rpeak     Rss      Tpeak    T50    ratio_Rpeak_Rss    Delta     L      peak_type
    _____    ______    _____    ___    _______________    _____    ____    _________

     NaN     8.5943     NaN     NaN          NaN           NaN     0.01      None   

PhaseTable = 1x8 table
    Rpeak      Rss      Tpeak       T50      ratio_Rpeak_Rss    Delta     L    peak_type
    ______    ______    ______    _______    _______________    ______    _    _________

    43.923    38.796    79.329    -53.719        1.1321         5.1266    1      High   

PhaseTable = 1x8 table
    Rpeak      Rss      Tpeak       T50      ratio_Rpeak_Rss      Delta       L    peak_type
    ______    ______    ______    _______    _______________    __________    _    _________

    51.168    51.168    500.33    0.36928           1           3.5032e-05    5      High   

PhaseTable = 1x8 table
    Rpeak      Rss      Tpeak       T50      ratio_Rpeak_Rss      Delta       L    peak_type
    ______    ______    ______    _______    _______________    __________    _    _________

    38.462    38.462    1999.9    -1.2205           1           4.5737e-05    1      High   

PhaseTable = 1x8 table
    Rpeak      Rss     Tpeak       T50      ratio_Rpeak_Rss    Delta    L    peak_type
    ______    _____    ______    _______    _______________    _____    _    _________

    54.771    52.24    2069.4    -54.078        1.0484         2.531    7      High   

PhaseTable = 1x8 table
    Rpeak      Rss      Tpeak       T50       ratio_Rpeak_Rss      Delta       L     peak_type
    ______    ______    ______    ________    _______________    __________    __    _________

    53.193    53.192    3001.8    -0.86721           1           3.2624e-05    10      High   

PhaseTable = 1x8 table
    Rpeak      Rss      Tpeak      T50      ratio_Rpeak_Rss      Delta        L      peak_type
    ______    ______    ______    ______    _______________    __________    ____    _________

    1.2225    1.2225    4001.9    1.7583           1           2.9831e-07    0.01      High   

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 - prev_end));
                        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 - prev_end));
                         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 - prev_end));
                         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, TauKb_ON, RLC)
    Kb = Kb_Min; % 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_Min;
            else
                % prev_end = PhaseTimes(j);
                if j < numel(RLC)
                    if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
                        Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
                    elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
                        Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
                    elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
                        Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
                    elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
                        Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
                    end
                else
                    if RLC(j-1) < RLC(j)
                        Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));

                    elseif RLC(j-1) > RLC(j)
                        Kb = Kb_Max - (Kb_Max - Kb_Min ) * exp(TauKb_ON * (T_end - T_start));
                    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

Furthermore, as @Torsten mentioned in his answer, you can instead use a 'for' loop function.

Hope this helps!

Answer 2

[Torsten]

Use

    PhaseTable = struct2table(PhaseResults);
    writetable(PhaseTable,'Phase.csv');

instead of

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

You only got one struct with name "PhaseResults".

And if you had several structs PhaseResults, you couldn't access them by "PhaseResults(i)", but by "PhaseResults{i}".

[Ehtisham]

PhaseTable = struct2table(PhaseResults);

writetable(PhaseTable,'Phase.csv'); Using this it again giving the csv file of the Phase1 just @Torsten

[Torsten]

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.01, 1, 5, 1, 7, 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, 3000, 4000,5000];% phase times (each row defines a phase)
timespan = [0, 5000];  % 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)

% Run the 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);

PhaseTable = 7x8 table
    Rpeak      Rss      Tpeak       T50       ratio_Rpeak_Rss      Delta        L      peak_type
    ______    ______    ______    ________    _______________    __________    ____    _________

       NaN    8.5943       NaN         NaN           NaN                NaN    0.01    {'None'} 
    43.923    38.796    79.329     -53.719        1.1321             5.1266       1    {'High'} 
    51.168    51.168    500.33     0.36928             1         3.5032e-05       5    {'High'} 
    38.462    38.462    1999.9     -1.2205             1         4.5737e-05       1    {'High'} 
    54.771     52.24    2069.4     -54.078        1.0484              2.531       7    {'High'} 
    53.193    53.192    3001.8    -0.86721             1         3.2624e-05      10    {'High'} 
    1.2225    1.2225    4001.9      1.7583             1         2.9831e-07    0.01    {'High'} 

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 maximum forward reaction rate
    Kf_LMax = Kf_Max * (RLC ./ (RLC + 1));
    
    for j = 1:size(PhaseTimes,1)
      % Calculate phase results
      for i = 1:size(PhaseTimes,2)-1
        PhaseResults{i} = calculate_phase(t, Active_Receptor_concentration, PhaseTimes(j,i:i+1), RLC(i));
      end
      PhaseResults = vertcat(PhaseResults{:});  % Concatenate results
    
      % Write Phase results to CSV
      PhaseTable = struct2table(PhaseResults)
      writetable(PhaseTable, ['Phase', num2str(j), '.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 - prev_end));
                        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 - prev_end));
                         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 - prev_end));
                         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, TauKb_ON, RLC)
    Kb = Kb_Min; % 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_Min;
            else
                % prev_end = PhaseTimes(j);
                if j < numel(RLC)
                    if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
                        Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
                    elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
                        Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
                    elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
                        Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
                    elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
                        Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
                    end
                else
                    if RLC(j-1) < RLC(j)
                        Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));

                    elseif RLC(j-1) > RLC(j)
                        Kb = Kb_Max - (Kb_Max - Kb_Min ) * exp(TauKb_ON * (T_end - T_start));
                    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

[Ehtisham]

One more question why after the 4000-5000 the Active receptor is not increasaing upto the value in the initial Phase at 0-10 @Torsten@Neelanshu