Why the if loop are getting the exact values of the Kf_LMax values not the approximated values in the different phase of the Relative ligand Concentration?

Why the if loop are getting the exact values of the Kf_LMax values not the approximated values in the different phase of the Relative ligand Concentration?

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 1, 5, 10, 0.1];  % RLC values
TauKf_ON = -0.1; % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF  
PhaseTimes = [0, 50, 100, 200, 400, 600, 1000]; % phase times (each row defines a phase)
% Generate time vector
t = 0:0.01:1000;
% Call the function to compute receptor concentrations and Kf
[ Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t);
% Plotting
figure;
% Plot Kf_LMax
plot(t, Kf_LMax, '-', 'LineWidth', 2 );
title('Kf_LMax');
xlabel('Time');
ylabel('Kf_LMax');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)
% Compute Kf_LMax values based on RLC
Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));

% Initialize output array
Kf_LMax = zeros(size(t));

% Number of phases
num_phases = numel(RLC);

% For each phase, calculate the corresponding Kf_LMax
for i = 1:num_phases
    % Get start and end times for the phase
    T_start = PhaseTimes(i);
    if i < num_phases
        T_end = PhaseTimes(i + 1);
    else
        T_end = inf;  % Last phase continues indefinitely
    end
    
    % Logical indices for time points in the current phase
    phase_idx = (t >= T_start) & (t < T_end);
    
    if i == 1
        % For the first phase, simply assign the maximum value
        Kf_LMax(phase_idx) = Kf_LMax_values(i);
    else
        % For later phases, handle the dynamic response
        prev_end = PhaseTimes(i - 1);
        if RLC(i - 1) < RLC(i)
            % Increasing RLC
            Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax_values(i - 1)) * exp(TauKf_ON * (t(phase_idx) - prev_end));
            
        else
            % Decreasing RLC
            Kf_LMax(phase_idx) = Kf_LMax_values(i) + (Kf_LMax_values(i) - Kf_LMax_values(i - 1)) * exp(TauKf_OFF * (t(phase_idx) - prev_end));
        end
    end
end
end

[dpb]

(Kf_LMax_values(i) - Kf_LMax_values(i - 1)) --> 0 because you set the fisrt values all to Kf_LMax_values

[Ehtisham]

@dpb i even did the modification but getting the same result

if i == 1
    % For the first phase, simply assign the maximum value
    Kf_LMax(phase_idx) = Kf_LMax_values(1);
    
else
    % For later phases, handle the dynamic response
    prev_end = PhaseTimes(i - 1);
    if RLC(i - 1) < RLC(i)
        % Increasing RLC
        Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax_values(i - 1)) * exp(TauKf_ON * (t(phase_idx) - prev_end));
        
    else
        % Decreasing RLC
        Kf_LMax(phase_idx) = Kf_LMax_values(i) + (Kf_LMax_values(i) - Kf_LMax_values(i - 1)) * exp(TauKf_OFF * (t(phase_idx) - prev_end));
    end
end
end

[dpb]

You're still multiplying the exponential portion by zero because the difference term is still there...see the Answer below.

Set a breakpoint at the line and hand calculate the terms and it will probably become clear...

[Voss]

"(Kf_LMax_values(i) - Kf_LMax_values(i - 1)) --> 0"

Not true.

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 1, 5, 10, 0.1];  % RLC values
TauKf_ON = -0.1; % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF  
PhaseTimes = [0, 50, 100, 200, 400, 600, 1000]; % phase times (each row defines a phase)
% Generate time vector
t = 0:0.01:1000;
% Call the function to compute receptor concentrations and Kf
[ Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t);

i=2: Kf_LMax_values(2) - Kf_LMax_values(1) = 24.242424
i=3: Kf_LMax_values(3) - Kf_LMax_values(2) = 16.666667
i=4: Kf_LMax_values(4) - Kf_LMax_values(3) = 33.333333
i=5: Kf_LMax_values(5) - Kf_LMax_values(4) = 7.575758
i=6: Kf_LMax_values(6) - Kf_LMax_values(5) = -81.818182

% Plotting
figure;
% Plot Kf_LMax
plot(t, Kf_LMax, '-', 'LineWidth', 2 );
title('Kf_LMax');
xlabel('Time');
ylabel('Kf_LMax');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)
% Compute Kf_LMax values based on RLC
Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));
% Initialize output array
Kf_LMax = zeros(size(t));
% Number of phases
num_phases = numel(RLC);
% For each phase, calculate the corresponding Kf_LMax
for i = 1:num_phases
    % Get start and end times for the phase
    T_start = PhaseTimes(i);
    if i < num_phases
        T_end = PhaseTimes(i + 1);
    else
        T_end = inf;  % Last phase continues indefinitely
    end
    
    % Logical indices for time points in the current phase
    phase_idx = (t >= T_start) & (t < T_end);
    
    if i == 1
        % For the first phase, simply assign the maximum value
        Kf_LMax(phase_idx) = Kf_LMax_values(i);
    else
        
        fprintf('i=%d: Kf_LMax_values(%d) - Kf_LMax_values(%d) = %f\n', ...
            i,i,i-1,Kf_LMax_values(i) - Kf_LMax_values(i - 1))
        
        % For later phases, handle the dynamic response
        prev_end = PhaseTimes(i - 1);
        if RLC(i - 1) < RLC(i)
            % Increasing RLC
            Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax_values(i - 1)) * exp(TauKf_ON * (t(phase_idx) - prev_end));
            
        else
            % Decreasing RLC
            Kf_LMax(phase_idx) = Kf_LMax_values(i) + (Kf_LMax_values(i) - Kf_LMax_values(i - 1)) * exp(TauKf_OFF * (t(phase_idx) - prev_end));
        end
    end
end
end

[Image Analyst]

To have underlines in the string be underlines and not cause a subscript, use the 'interpreter' 'none' option

title('Kf_LMax', 'Interpreter', 'none');
xlabel('Time', 'Interpreter', 'none');
ylabel('Kf_LMax', 'Interpreter', 'none');

Answer 1

[Voss]

prev_end = PhaseTimes(i - 1);

That's the start of the previous phase, not the end of the previous phase.

Using the previous phase's start time causes the exponentials to appear to be approximately flat by the time the current phase starts, which is why the plot looks like a bunch of horizontal lines.

My guess is that you want to see exponential ramp up in each phase (except the first phase), and I guess for that you should use the end of the previous phase (= start of the current phase = PhaseTimes(i) = T_start).

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 1, 5, 10, 0.1];  % RLC values
TauKf_ON = -0.1; % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF  
PhaseTimes = [0, 50, 100, 200, 400, 600, 1000]; % phase times (each row defines a phase)
% Generate time vector
t = 0:0.01:1000;
% Call the function to compute receptor concentrations and Kf
[ Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t);
% Plotting
figure;
% Plot Kf_LMax
plot(t, Kf_LMax, '-', 'LineWidth', 2 );
title('Kf_LMax');
xlabel('Time');
ylabel('Kf_LMax');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)
% Compute Kf_LMax values based on RLC
Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));
% Initialize output array
Kf_LMax = zeros(size(t));
% Number of phases
num_phases = numel(RLC);
% For each phase, calculate the corresponding Kf_LMax
for i = 1:num_phases
    % Get start and end times for the phase
    T_start = PhaseTimes(i);
    if i < num_phases
        T_end = PhaseTimes(i + 1);
    else
        T_end = inf;  % Last phase continues indefinitely
    end
    
    % Logical indices for time points in the current phase
    phase_idx = (t >= T_start) & (t < T_end);
    
    if i == 1
        % For the first phase, simply assign the maximum value
        Kf_LMax(phase_idx) = Kf_LMax_values(i);
    else
        
        % For later phases, handle the dynamic response
        if RLC(i - 1) < RLC(i)
            % Increasing RLC
            Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax_values(i - 1)) * exp(TauKf_ON * (t(phase_idx) - T_start));
            
        else
            % Decreasing RLC
            Kf_LMax(phase_idx) = Kf_LMax_values(i) + (Kf_LMax_values(i) - Kf_LMax_values(i - 1)) * exp(TauKf_OFF * (t(phase_idx) - T_start));
        end
    end
end
end

[Ehtisham]

@Voss @Image AnalystI used that code that written above but they giving the shoot up and even though it is a exponentinal function see the picture. How i can get rid of these shoot up

[Image Analyst]

I don't know. It's not my code. What/where would you like the signal to be right after it jumps? Lowered so that it's almost continuous with the prior signal? If so, why? What would happen to the very peak? If not lowered, then where would the new signal lie?

[Voss]

@Ehtisham: The code in my answer generates the plot in my answer. If you got a different plot, then you used a different code. Please share the code you used so that people can reproduce your plot, and then please answer @Image Analyst's questions about how/why you want the plot to be different.

[Ehtisham]

@Voss I just changes the Tau values to see the response of the curve but its shoot up and not starting the new phase in the end of the last phase. @Image Analyst it is a exponential function so it should start the new phase in the end of the last phase.

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 1, 2, 5, 0.1];  % RLC values
TauKf_ON = -0.01; % TauKf_ON
TauKf_OFF = -0.01;  % TauKf_OFF  
PhaseTimes = [0, 100, 200, 300, 400, 500, 600]; % phase times (each row defines a phase)
% Generate time vector
t = 0:0.1:600;
% Call the function to compute receptor concentrations and Kf
[ Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t);
% Plotting
figure;
% Plot Kf_LMax
plot(t, Kf_LMax, 'b', 'LineWidth', 2 );
title('Kf_LMax');
xlabel('Time');
ylabel('Kf_LMax');
grid off;

function [Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)
% Compute Kf_LMax values based on RLC
Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));  % Normalized Kf_LMax values
% Initialize output array
Kf_LMax = zeros(size(t));  % Array to store calculated values
% Number of phases
num_phases = numel(RLC);
% Loop over each phase
for i = 1:num_phases
    % Get start and end times for each phase
    T_start = PhaseTimes(i);
    if i < num_phases
        T_end = PhaseTimes(i + 1);
    else
        T_end = inf;  % Last phase continues indefinitely
    end
    % Identify the indices for the current phase time range
    phase_idx = (t >= T_start) & (t < T_end);
    if i == 1
        % In the first phase, assign Kf_LMax as the calculated value
        Kf_LMax(phase_idx) = Kf_LMax_values(i);
    else
        
        if RLC(i - 1) < RLC(i)
            % If RLC is increasing, apply TauKf_ON for a smooth transition
            Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax_values(i - 1)) * exp(TauKf_ON * (t(phase_idx) - T_start));
        else
            % If RLC is decreasing, apply TauKf_OFF for a smooth transition
            Kf_LMax(phase_idx) = Kf_LMax_values(i) + (Kf_LMax_values(i - 1) - Kf_LMax_values(i)) * exp(TauKf_OFF * (t(phase_idx) - T_start));
        end
    end
end
end

[Voss]

When tau=-0.1, each exponential rise (or fall) reaches a level that is approximately its steady state level (i.e., the level at t=infinity), so using that level or using the actual steady-state value (Kf_LMax_values(i-1)) makes no visual difference to the plot.

However, when tau=-0.01, the exponentials don't have enough time to reach near their steady-state levels, so you get those jumps in the plot, where the next phase starts at the prescribed level but the previous phase didn't reach it.

To use the actual last level reached of the previous phase as the starting level of the each phase, replace Kf_LMax_values(i-1) with Kf_LMax(last_phase_idx), where last_phase_idx is the index in t of the final point of the previous phase.

Making that change shows that the plot when tau=-0.1 is the same as before:

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 1, 2, 5, 0.1];  % RLC values
TauKf_ON = -0.1; % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF  
PhaseTimes = [0, 100, 200, 300, 400, 500, 600]; % phase times (each row defines a phase)
% Generate time vector
t = 0:0.1:600;
% Call the function to compute receptor concentrations and Kf
[ Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t);
% Plotting
figure;
% Plot Kf_LMax
plot(t, Kf_LMax, 'b', 'LineWidth', 2 );
title('Kf\_LMax');
xlabel('Time');
ylabel('Kf\_LMax');
grid off;

But the plot with tau=-0.01 is now without the jumps:

TauKf_ON = -0.01; % TauKf_ON
TauKf_OFF = -0.01;  % TauKf_OFF  
% Call the function to compute receptor concentrations and Kf
[ Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t);
% Plotting
figure;
% Plot Kf_LMax
plot(t, Kf_LMax, 'b', 'LineWidth', 2 );
title('Kf\_LMax');
xlabel('Time');
ylabel('Kf\_LMax');
grid off;

function [Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)
% Compute Kf_LMax values based on RLC
Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));  % Normalized Kf_LMax values
% Initialize output array
Kf_LMax = zeros(size(t));  % Array to store calculated values
% Number of phases
num_phases = numel(RLC);
% Loop over each phase
for i = 1:num_phases
    % Get start and end times for each phase
    T_start = PhaseTimes(i);
    if i < num_phases
        T_end = PhaseTimes(i + 1);
    else
        T_end = inf;  % Last phase continues indefinitely
    end
    % Identify the indices for the current phase time range
    phase_idx = (t >= T_start) & (t < T_end);
    if i == 1
        % In the first phase, assign Kf_LMax as the calculated value
        Kf_LMax(phase_idx) = Kf_LMax_values(i);
    else
        if RLC(i - 1) < RLC(i)
            % If RLC is increasing, apply TauKf_ON for a smooth transition
            Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) * exp(TauKf_ON * (t(phase_idx) - T_start));
        else
            % If RLC is decreasing, apply TauKf_OFF for a smooth transition
            Kf_LMax(phase_idx) = Kf_LMax_values(i) + (Kf_LMax(last_phase_idx) - Kf_LMax_values(i)) * exp(TauKf_OFF * (t(phase_idx) - T_start));
        end
    end
    last_phase_idx = find(phase_idx,1,'last');
end
end

[Voss]

[Replying to a deleted comment by @Ehtisham]

@Ehtisham: The problem seems to be that you changed some of the variable names (Kf_LMax became Kf_L, and Kf_LMax_values became Kf_LMax) and made a mistake when you did so: you have Kf_LMax(last_phase_idx) instead of Kf_L(last_phase_idx). [The idea is to take the value of the curve you're calculating, at the end of the previous phase, as the reference level, so that's Kf_L at last_phase_idx.]

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 1, 2, 5, 0.1];  % RLC values
TauKf_ON = -0.01; % TauKf_ON
TauKf_OFF = -0.01;  % TauKf_OFF  
PhaseTimes = [0, 100, 200, 300, 400, 500, 600]; % phase times (each row defines a phase)
% Generate time vector
t = 0:0.1:600;
% Call the function to compute receptor concentrations and Kf
[ Kf_LMax] = calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);
% Plotting
figure;
% Plot Kf_LMax
plot(t, Kf_LMax, 'b', 'LineWidth', 2 );
title('Kf\_LMax');
xlabel('Time');
ylabel('Kf\_LMax');
grid off;

function Kf_L = calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)
% Compute Kf_LMax values based on RLC
Kf_LMax = Kf_Max * (RLC ./ (1 + RLC));  % Normalized Kf_LMax values
% Initialize output array
Kf_L = zeros(size(t));  % Array to store calculated values
% Number of phases
num_phases = numel(RLC);
% Loop over each phase
for i = 1:num_phases
    % Get start and end times for each phase
    T_start = PhaseTimes(i);
    if i < num_phases
        T_end = PhaseTimes(i + 1);
    else
        T_end = inf;  % Last phase continues indefinitely
    end
    % Identify the indices for the current phase time range
    phase_idx = (t >= T_start) & (t < T_end);
    if i == 1
        % In the first phase, assign Kf_LMax as the calculated value
        Kf_L(phase_idx) = Kf_LMax(i);
    else
        if RLC(i - 1) < RLC(i)
            % If RLC is increasing, apply TauKf_ON for a smooth transition
            Kf_L(phase_idx) = Kf_LMax(i) - (Kf_LMax(i) - Kf_L(last_phase_idx)) * exp(TauKf_ON * (t(phase_idx) - T_start));
        else
            % If RLC is decreasing, apply TauKf_OFF for a smooth transition
            Kf_L(phase_idx) = Kf_LMax(i) + (Kf_L(last_phase_idx) - Kf_LMax(i)) * exp(TauKf_OFF * (t(phase_idx) - T_start));
        end
    end
    last_phase_idx = find(phase_idx,1,'last');
end 
end

[Sam Chak]

Hi @Ehtisham,

To be honest, how confident are you in accurately modeling the continuous-time dynamics of the Relative Ligand Concentration using multi-level or nested If–Else statements? Are you aware of the potential pitfalls?

Additionally, if you plan to publish the results in reputable journals, will you describe the changes in ligand molecules over time using differential equations or nested If–Else statements? These questions are important because the editor or reviewers may raise concerns about the validity of the methodology and results.

[Voss]

@Ehtisham: If my answer satisfactorily answers your question, please Accept it. Thanks!

[Torsten]

Where and how do you call "SimfileNPhase" ?

With the included code, no such error happens.

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 1, 5, 10, 0.1];  % RLC values
TauKf_ON = -0.1; % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF  
PhaseTimes = [0, 50, 100, 200, 400, 600, 1000]; % phase times (each row defines a phase)
% Generate time vector
t = 0:0.01:1000;
% Call the function to compute receptor concentrations and Kf
[ Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t);
% Plotting
figure;
% Plot Kf_LMax
plot(t, Kf_LMax, '-', 'LineWidth', 2 );
title('Kf_LMax');
xlabel('Time');
ylabel('Kf_LMax');
grid on;

function  [t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase (Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...
    Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, PhaseTimes,timespan,offsets,Non_Active_Receptor_concentration, Active_Receptor_concentration)

    % Define functions for forward and backward reaction rates
    Kf_L = @(t) Kf_Cal (t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);
    Kb = @(t) Kb_Cal(Kb_Max, Kb_Min, RLC, TauKb_ON, TauKb_OFF, PhaseTimes, t, offsets);

    % 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.01, 'RelTol', 1e-4, 'AbsTol', 1e-6);

    % 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 / (1 + RLC));

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Nested function for calculating Kf_L based on time
function [Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)
% Compute Kf_LMax values based on RLC
    % Compute Kf_LMax values based on RLC
    Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));  % Normalized Kf_LMax values
    % Initialize output array
    Kf_LMax = zeros(size(t));  % Array to store calculated values
    % Number of phases
    num_phases = numel(RLC);
    % Loop over each phase
    for i = 1:num_phases
        % Get start time for each phase
        T_start = PhaseTimes(i);
        % Get end time for the current phase, or set to infinity if it's the last phase
        if i < num_phases
            T_end = PhaseTimes(i + 1);
        else
            T_end = inf;  % Last phase continues indefinitely
        end
       
        % Identify the indices for the current phase time range
        phase_idx = (t >= T_start) & (t < T_end);
        
        if i == 1
            % In the first phase, assign Kf_LMax as the calculated value
            Kf_LMax(phase_idx) = Kf_LMax_values(i);
        else
            if RLC(i - 1) < RLC(i)
                % If RLC is increasing, apply TauKf_ON for a smooth transition
                Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) * exp(TauKf_ON * (t(phase_idx) - T_start));
            else
                % If RLC is decreasing, apply TauKf_OFF for a smooth transition
                Kf_LMax(phase_idx) = Kf_LMax_values(i) + (Kf_LMax(last_phase_idx) - Kf_LMax_values(i)) * exp(TauKf_OFF * (t(phase_idx) - T_start));
            end
        end
        % Save the last valid index for smooth transitions between phases
        last_phase_idx = find(phase_idx, 1, 'last');
    end
end

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kb = Kb_Cal(Kb_Max, Kb_Min, RLC, TauKb_ON, TauKb_OFF, PhaseTimes, t, offsets)
    % Determine which phase each time step falls into
    num_phases = numel(RLC);
    phase_idx = zeros(size(t));
    
    for j = 1:num_phases
        T_start = PhaseTimes(j);
        if j < num_phases
            T_end = PhaseTimes(j + 1);
        else
            T_end = inf;
        end
        phase_idx(t >= T_start & t < T_end) = j;
    end

    % Ensure offsets are provided for each phase
    if numel(offsets) ~= num_phases
        error('Number of offsets must match the number of phases in RLC.');
    end

    % Calculate Kb for each time step based on the phase
    Kb = zeros(size(t));
    for i = 1:length(t)
        j = phase_idx(i);  % Determine current phase
        current_offset = offsets(j);  % Use the offset for the current phase
        prev_end = PhaseTimes(max(j - 1, 1));  % Get the previous phase end time
        
        if j == 1
            % For the first phase, start at Kb_Max
            Kb(i) = Kb_Max;
        else
            % Calculate Kb depending on whether RLC is increasing or decreasing
            if RLC(j - 1) < RLC(j)
                % If RLC increases, use TauKb_ON
                Kb(i) = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t(i) - prev_end)) + current_offset;
            else
                % If RLC decreases, use TauKb_OFF
                Kb(i) = Kb_Min + (Kb_Max - Kb_Min) * exp(TauKb_OFF * (t(i) - prev_end)) + current_offset;
            % Ensure Kb does not exceed Kb_Max
               Kb(end)= Kb_Max;
            end
            
        end
       
    end
  
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Ensure the ode_LR function is properly set to work with 2 elements in y:
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(Phase_time, Phase_concentration)
    % Compute the derivative 
    dt = diff(Phase_time);
    dy = diff(Phase_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(Phase_concentration(max_indices));
                % Find the lowest minimum
                [Rpeak_min, minIdx] = min(Phase_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 = Phase_time(max_indices(maxIdx));
                    peak_type = 'High';
                else
                    Rpeak = Rpeak_min;
                    Tpeak = Phase_time(min_indices(minIdx));
                    peak_type = 'Low';
                end
            % If only maxima exist
            elseif ~isempty(max_indices)
                [Rpeak, maxIdx] = max(Phase_concentration(max_indices));
                Tpeak = Phase_time(max_indices(maxIdx));
                peak_type = 'High';
            % If only minima exist
            elseif ~isempty(min_indices)
                [Rpeak, minIdx] = min(Phase_concentration(min_indices));
                Tpeak = Phase_time(min_indices(minIdx));
                peak_type = 'Low';
            end
        end
    end
end

[Torsten]

The list of input arguments to Kf_Cal isn't consistent.

Kf_L = @(t) Kf_Cal (t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);

function [Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)

[Ehtisham]

@Torsten I fix that issuse but now am getting this issue

Unable to perform assignment because the left and right sides have a different number of elements.

Error in SimfileNPhase>Kf_Cal (line 96)

Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) * exp(TauKf_ON * (t(phase_idx) - T_start));

Error in SimfileNPhase>@(t)Kf_Cal(t,PhaseTimes,Kf_Max,RLC,TauKf_ON,TauKf_OFF) (line 5)

Kf_L = @(t) Kf_Cal (t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);

Error in SimfileNPhase>ode_LR (line 245)

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

Error in SimfileNPhase>@(t,y)ode_LR(t,y,Kf_L,Kb) (line 19)

[t, y] = ode45(@(t, y) ode_LR(t, y, Kf_L, Kb), timespan, initial_conditions, options);

Error in ode45 (line 300)

f4 = ode(t4, y4);

Error in SimfileNPhase (line 19)

[t, y] = ode45(@(t, y) ode_LR(t, y, Kf_L, Kb), timespan, initial_conditions, options);

Error in SimfileNPhaseRun (line 21)

[t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase( ...

>>

[Torsten]

It seems you deleted the script part to reproduce the error.

[Ehtisham]

The script is same . here the attach file @Torsten

[Torsten]

I miss the calling (script) part for your functions.

The part that replaces

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 1, 5, 10, 0.1];  % RLC values
TauKf_ON = -0.1; % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF  
PhaseTimes = [0, 50, 100, 200, 400, 600, 1000]; % phase times (each row defines a phase)
% Generate time vector
t = 0:0.01:1000;
% Call the function to compute receptor concentrations and Kf
[ Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t);
% Plotting
figure;
% Plot Kf_LMax
plot(t, Kf_LMax, '-', 'LineWidth', 2 );
title('Kf_LMax');
xlabel('Time');
ylabel('Kf_LMax');
grid on;

[Ehtisham]

@Torsten % Define parameters

Kf_Max = 100; % maximum forward rate

Kb_Max = 80; % maximum backward rate

Kb_Min = 10;

RLC = [0.1, 0.5, 1, 10, 5, 0.1]; % 4 RLC values for 4 phases

offsets = [10, 20, 30, 40, 50, 60];

TauKf_ON = -0.1; % TauKf_ON

TauKf_OFF = -0.1; % TauKf_OFF

TauKb_ON = -0.1; % TauKb_ON

TauKb_OFF = -0.1; % TauKb_OFF

% PhaseTimes should have one more element than RLC

PhaseTimes = [0, 500, 1000, 2000, 2500, 3000];

timespan = 0:0.01:3000; % Adjusted timespan to match the extended PhaseTimes

% Example initial conditions for non-active and active receptor concentrations

Non_Active_Receptor_concentration = 100; % Initial non-active receptor concentration

Active_Receptor_concentration = 0; % Initial active receptor concentration

% Call the simulation function

[t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase( ...

Kf_Max, RLC, TauKf_ON, TauKf_OFF, Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, ...

PhaseTimes, timespan, offsets, Non_Active_Receptor_concentration, Active_Receptor_concentration);

[Torsten]

As you can see from the code, last_phase_idx becomes empty. In this case,

Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) .* exp(TauKf_ON * (t(phase_idx) - T_start))

becomes empty, and the assignment

Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) .* exp(TauKf_ON * (t(phase_idx) - T_start))

no longer works.

Further, phase_idx is always either logical 0 or logical 1. Thus I don't understand the assignment in total.

% Define parameters
Kf_Max = 100;  % maximum forward rate
Kb_Max = 80;  % maximum backward rate 
Kb_Min = 10;
RLC = [0.1, 0.5, 1, 10, 5, 0.1];  % 4 RLC values for 4 phases
offsets = [10, 20, 30, 40, 50, 60];
TauKf_ON = -0.1;  % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF
TauKb_ON = -0.1;  % TauKb_ON
TauKb_OFF = -0.1;  % TauKb_OFF
% PhaseTimes should have one more element than RLC
PhaseTimes = [0, 500, 1000, 2000, 2500, 3000];  
timespan = 0:0.01:3000;  % Adjusted timespan to match the extended PhaseTimes
% Example initial conditions for non-active and active receptor concentrations
Non_Active_Receptor_concentration = 100;  % Initial non-active receptor concentration
Active_Receptor_concentration = 0;  % Initial active receptor concentration
% Call the simulation function
[t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase( ...
    Kf_Max, RLC, TauKf_ON, TauKf_OFF, Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, ...
    PhaseTimes, timespan, offsets, Non_Active_Receptor_concentration, Active_Receptor_concentration);

phase_idx = logical
   1

last_phase_idx =

     []

s1 = 1×2

1 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

s2 = 1×2

0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Unable to perform assignment because the left and right sides have a different number of elements.

Error in solution>Kf_Cal (line 124)
Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) .* exp(TauKf_ON * (t(phase_idx) - T_start));

Error in solution>@(t)Kf_Cal(t,PhaseTimes,Kf_Max,RLC,TauKf_ON,TauKf_OFF) (line 25)
Kf_L = @(t) Kf_Cal (t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);

Error in solution>ode_LR (line 238)
dNonActiveReceptor_dt = -Kf_L(t) * Non_Active_Receptor_concentration + Kb(t) * Active_Receptor_concentration;

Error in solution>@(t,y)ode_LR(t,y,Kf_L,Kb) (line 39)
[t, y] = ode45(@(t, y) ode_LR(t, y, Kf_L, Kb), timespan, initial_conditions, options);

Error in ode45 (line 293)
f4 = ode(t4, y4, odeArgs{:});

Error in solution>SimfileNPhase (line 39)
[t, y] = ode45(@(t, y) ode_LR(t, y, Kf_L, Kb), timespan, initial_conditions, options);

function  [t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase (Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...
    Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, PhaseTimes,timespan,offsets,Non_Active_Receptor_concentration, Active_Receptor_concentration)

    % Define functions for forward and backward reaction rates
    Kf_L = @(t) Kf_Cal (t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);
    Kb = @(t) Kb_Cal(Kb_Max, Kb_Min, RLC, TauKb_ON, TauKb_OFF, PhaseTimes, t, offsets);

    % 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.01, 'RelTol', 1e-4, 'AbsTol', 1e-6);

    % 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 / (1 + RLC));

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Nested function for calculating Kf_L based on time
function [Kf_LMax] = Kf_Cal(t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)
   % Compute Kf_LMax values based on RLC
    Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));  % Normalized Kf_LMax values
    % Initialize output array
    Kf_LMax = zeros(size(t));  % Array to store calculated values
    % Number of phases
    num_phases = numel(RLC);
    % Loop over each phase
    for i = 1:num_phases
        % Get start time for each phase
        T_start = PhaseTimes(i);
        % Get end time for the current phase, or set to infinity if it's the last phase
        if i < num_phases
            T_end = PhaseTimes(i + 1);
        else
            T_end = inf;  % Last phase continues indefinitely
        end

        % Identify the indices for the current phase time range
        phase_idx = (t >= T_start) & (t < T_end);

        if i == 1
            % In the first phase, assign Kf_LMax as the calculated value
            Kf_LMax(phase_idx) = Kf_LMax_values(i);
        else
            if RLC(i - 1) < RLC(i)
                % If RLC is increasing, apply TauKf_ON for a smooth transition
                s1 = size(Kf_LMax(phase_idx));
                s2 = size(Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) .* exp(TauKf_ON * (t(phase_idx) - T_start)));
                if s1 ~= s2
                  phase_idx
                  last_phase_idx
                  s1
                  s2
                end
                Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) .* exp(TauKf_ON * (t(phase_idx) - T_start));
            else
                % If RLC is decreasing, apply TauKf_OFF for a smooth transition
                Kf_LMax(phase_idx) = Kf_LMax_values(i) + (Kf_LMax(last_phase_idx) - Kf_LMax_values(i)) .* exp(TauKf_OFF * (t(phase_idx) - T_start));
            end
        end
        % % Save the last valid index for smooth transitions between phases
         last_phase_idx = find(phase_idx, 1, 'last');
    end
end
%     % Compute Kf_LMax values based on RLC
%     Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));  % Normalized Kf_LMax values
% 
%     % Initialize output array
%     Kf_LMax = zeros(size(t));  % Array to store calculated values
% 
%     % Number of phases
%     num_phases = numel(RLC);
% for i = 1:num_phases
%     T_start = PhaseTimes(i);
% 
%     if i < num_phases
%         T_end = PhaseTimes(i + 1);
%     else
%         T_end = inf;
%     end
% 
%     phase_idx = (t >= T_start) & (t < T_end);
% 
%     % Initialize last_phase_idx
%     if any(phase_idx)
%         last_phase_idx = find(phase_idx, 1, 'last');
%     else
%         last_phase_idx = NaN; 
%     end
% 
%     if i == 1
%         Kf_LMax(phase_idx) = Kf_LMax_values(i);
%     elseif i == num_phases
%         Kf_LMax(phase_idx) = Kf_LMax_values(i);
%     else
%         if isnan(last_phase_idx)
%             Kf_LMax(phase_idx) = Kf_LMax_values(i); % or handle accordingly
%         else
%             if RLC(i - 1) < RLC(i)
%                 Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) * exp(TauKf_ON * (t(phase_idx) - T_start));
%             else
%                 Kf_LMax(phase_idx) = Kf_LMax_values(i) + (Kf_LMax(last_phase_idx) - Kf_LMax_values(i)) * exp(TauKf_OFF * (t(phase_idx) - T_start));
%             end
%         end
%     end
% end
% 
% 
% end


 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kb = Kb_Cal(Kb_Max, Kb_Min, RLC, TauKb_ON, TauKb_OFF, PhaseTimes, t, offsets)
    % Determine which phase each time step falls into
    num_phases = numel(RLC);
    phase_idx = zeros(size(t));
    
    for j = 1:num_phases
        T_start = PhaseTimes(j);
        if j < num_phases
            T_end = PhaseTimes(j + 1);
        else
            T_end = inf;
        end
        phase_idx(t >= T_start & t < T_end) = j;
    end

    % Ensure offsets are provided for each phase
    if numel(offsets) ~= num_phases
        error('Number of offsets must match the number of phases in RLC.');
    end

    % Calculate Kb for each time step based on the phase
    Kb = zeros(size(t));
    for i = 1:length(t)
        j = phase_idx(i);  % Determine current phase
        current_offset = offsets(j);  % Use the offset for the current phase
        prev_end = PhaseTimes(max(j - 1, 1));  % Get the previous phase end time
        
        if j == 1
            % For the first phase, start at Kb_Max
            Kb(i) = Kb_Max;
        else
            % Calculate Kb depending on whether RLC is increasing or decreasing
            if RLC(j - 1) < RLC(j)
                % If RLC increases, use TauKb_ON
                Kb(i) = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t(i) - prev_end)) + current_offset;
            else
                % If RLC decreases, use TauKb_OFF
                Kb(i) = Kb_Min + (Kb_Max - Kb_Min) * exp(TauKb_OFF * (t(i) - prev_end)) + current_offset;
            % Ensure Kb does not exceed Kb_Max
               Kb(end)= Kb_Max;
            end
            
        end
       
    end
  
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Ensure the ode_LR function is properly set to work with 2 elements in y:
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(Phase_time, Phase_concentration)
    % Compute the derivative 
    dt = diff(Phase_time);
    dy = diff(Phase_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(Phase_concentration(max_indices));
                % Find the lowest minimum
                [Rpeak_min, minIdx] = min(Phase_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 = Phase_time(max_indices(maxIdx));
                    peak_type = 'High';
                else
                    Rpeak = Rpeak_min;
                    Tpeak = Phase_time(min_indices(minIdx));
                    peak_type = 'Low';
                end
            % If only maxima exist
            elseif ~isempty(max_indices)
                [Rpeak, maxIdx] = max(Phase_concentration(max_indices));
                Tpeak = Phase_time(max_indices(maxIdx));
                peak_type = 'High';
            % If only minima exist
            elseif ~isempty(min_indices)
                [Rpeak, minIdx] = min(Phase_concentration(min_indices));
                Tpeak = Phase_time(min_indices(minIdx));
                peak_type = 'Low';
            end
        end
    end
end

[Ehtisham]

@Voss Need your guidance in this matter. Why i am getting this error

Unable to perform assignment because the left and right sides have a different number of elements.

Error in SimfileNPhase>Kf_Cal (line 97) Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) * exp(TauKf_ON * (t(phase_idx) - T_start));

[Torsten]

But I gave you the answer: Because "last_phase_idx" becomes empty in the course of your calculations.

[Walter Roberson]

@Ehtisham

@Torsten explains the reasons for that error: last_phase_idx becomes empty. In any case, last_phase_idx is logical and the assignment makes no sense when last_phase_idx is false.

[Ehtisham]

@Torsten @Walter Roberson How i can get rid of this error

[Walter Roberson]

After

phase_idx = (t >= T_start) & (t < T_end);

insert

if ~any(phase_idx); continue; end

[Ehtisham]

Unrecognized function or variable 'last_phase_idx'.

Error in SimfileNPhase>Kf_Cal (line 97)

Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) * exp(TauKf_ON * (t(phase_idx) - T_start));

Error in SimfileNPhase>@(t)Kf_Cal(t,PhaseTimes,Kf_Max,RLC,TauKf_ON,TauKf_OFF) (line 5)

Kf_L = @(t) Kf_Cal (t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);

Now getting this error @Walter Roberson

function [Kf_LMax] = Kf_Cal(t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)
% Compute Kf_LMax values based on RLC
%     Compute Kf_LMax values based on RLC
Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));  % Normalized Kf_LMax values
% Initialize output array
Kf_LMax = zeros(size(t));  % Array to store calculated values
%     Number of phases
num_phases = numel(RLC);
%     Loop over each phase
for i = 1:num_phases
    % Get start time for each phase
    T_start = PhaseTimes(i);
    %         Get end time for the current phase, or set to infinity if it's the last phase
    if i < num_phases
        T_end = PhaseTimes(i + 1);
    else
        T_end = inf;  % Last phase continues indefinitely
    end
    % 
    %         Identify the indices for the current phase time range
    phase_idx = (t >= T_start) & (t < T_end);
    if ~any(phase_idx); continue; end
    if i == 1
        %             In the first phase, assign Kf_LMax as the calculated value
        Kf_LMax(phase_idx) = Kf_LMax_values(i);
    else
        
        if RLC(i - 1) < RLC(i)
            % If RLC is increasing, apply TauKf_ON for a smooth transition
            Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) * exp(TauKf_ON * (t(phase_idx) - T_start));
        else
            % If RLC is decreasing, apply TauKf_OFF for a smooth transition
            Kf_LMax(phase_idx) = Kf_LMax_values(i) + (Kf_LMax(last_phase_idx) - Kf_LMax_values(i)) * exp(TauKf_OFF * (t(phase_idx) - T_start));
        end
    end
    % Save the last valid index for smooth transitions between phases
    last_phase_idx = find(phase_idx, 1, 'last');
end

end

[Walter Roberson]

Suppose when i is 1, that ~any(phase_idx) is true. You will continue to the for i = 2 iteration.

Suppose (for the sake of argument) that ~any(phase_idx) is false when i is 2. Then if i == 1 will be false because i is 2. You proceed on to if RLC(i - 1) < RLC(i) . Now, no matter which of the two branches is chosen, the code needs last_phase_idx to be defined. But last_phase_idx was not defined, because you skipped the assignment to last_phase_idx due to the continue statement.

The test if i == 1 is not correct: what you need to know is not whether i is 1, but rather whether this is the first iteration in which phase_idx was usable.

[Ehtisham]

@Walter Roberson How i can solve this error

[Torsten]

We can only list the problems with your code, but cannot give advice how to solve the errors because we don't understand what you are trying to do.

[Walter Roberson]

@Ehtisham

You initialize a variable

last_phase_valid = false;

You change

    if i == 1

to

    if ~last_phase_valid

After

        last_phase_idx = find(phase_idx, 1, 'last');

you add

        last_phase_valid = true;

[Ehtisham]

function [Kf_LMax] = Kf_Cal(t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)
% Compute Kf_LMax values based on RLC
Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));  % Normalized Kf_LMax values
% Initialize output array
Kf_LMax = zeros(size(t));  % Array to store calculated values
% Number of phases
num_phases = numel(RLC);
% Initialize the last phase index and validity flag for transitions
last_phase_idx = 1;  % Initial index for the start of the first phase
last_phase_valid = false;  % Validity flag for transitions
% Loop over each phase
for i = 1:num_phases
    % Get start time for each phase
    T_start = PhaseTimes(i);
    % Get end time for the current phase, or set to infinity if it's the last phase
    if i < num_phases
        T_end = PhaseTimes(i + 1);
    else
        T_end = inf;  % Last phase continues indefinitely
    end
    % Identify the indices for the current phase time range
    phase_idx = (t >= T_start) & (t < T_end);
    % Transition handling for first and subsequent phases
    if ~last_phase_valid
        % In the first phase, assign Kf_LMax as the calculated value directly
        Kf_LMax(phase_idx) = Kf_LMax_values(i);
    else
        % Smooth transition based on the change in RLC value
        if RLC(i - 1) < RLC(i)
            % If RLC is increasing, apply TauKf_ON for a smooth transition
            Kf_LMax(phase_idx) = Kf_LMax_values(i) - ...
                (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) .* exp(TauKf_ON * (t(phase_idx) - T_start));
        else
            % If RLC is decreasing, apply TauKf_OFF for a smooth transition
            Kf_LMax(phase_idx) = Kf_LMax_values(i) + ...
                (Kf_LMax(last_phase_idx) - Kf_LMax_values(i)) .* exp(TauKf_OFF * (t(phase_idx) - T_start));
        end
    end
    % Save the last valid index for smooth transitions between phases
    last_phase_idx = find(phase_idx, 1, 'last');
    last_phase_valid = true;
end
end 

Unable to perform assignment because the left and right sides have a different number of elements.

Error in SimfileNPhase>Kf_Cal (line 119)

Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax(last_phase_idx)) .* exp(TauKf_ON .* (t(phase_idx) - T_start));

@Walter Roberson still getting the error

[Walter Roberson]

After

    phase_idx = (t >= T_start) & (t < T_end);

you still need

    if ~any(phase_idx); continue; end

[Ehtisham]

@Walter Roberson

Why i do not get the response of the Kf when i change the Tau values to 0.1 or 0.01 like picture below in my

and

code after the changes that you mention

[Walter Roberson]

I am going to need full current source from you to have a chance of answering this.

[Ehtisham]

@Walter Roberson

In Sim File my Kf code was giving the fixed values for app phases( like for RLC = 0.1 gives Kf values at 33.33 ) even when i change the Tau values so then but the code i put a question here and Torsten said your need to chnage the Kf_LMax_values(i-1) to Kf_LMax( last_phase_idx) that code work i give the graph that given above but when i put the Kf function in the Simfile function (main function ) its giving me assgiment error. after solving the assigmnet error again the Sime file (Kf function ) giving me fixed values like before.

[Sam Chak]

@Ehtisham,

You may have forgotten my previous comment, but this issue has persisted for more than three months. Please understand that eliminating the errors in the MATLAB code does not necessarily mean that your solution to the actual dynamical problem is correct. To address the problem effectively, you need to fully understand the mathematics behind it so that you can control the behavior of the system.

The exponential behaviors suggest that the system can be described by a 1st-order linear differential equation , in which the time constant of the system (possibly related to 'TauKf') influences the rate of exponential decay. As a researcher (not merely a programmer), you should have a single master mathematical model in the code, where you can flexibly set the value of the time constant.

Having too many nested 'If–Else' structures may hinder your ability to troubleshoot the dynamics easily, as each conditional statement can disrupt your visualization of the continuous time flow of the system. My suggestion is to apply a 'Divide and Conquer' approach by working through the time segments separately. If the exponential behaviors are correct, you can then combine them.

First and foremost, present the mathematical model, or, better yet, ask a Question. Do not simply accept a solution if it has not truly addressed your underlying problem.

Answer 2

[dpb]

Please format your code with the CODE button (or select and Ctrl-E)...

You're still multiplying the exponential portion by zero because the difference term is still there...probably what you're looking for is more like

Kf_LMax(phase_idx) = Kf_LMax_values(i-1)-Kf_LMax_values(i-1))*exp(TauKf_ON * (t(phase_idx) - prev_end));

which could be rewritten as

Kf_LMax(phase_idx)=Kf_LMax_values(i-1)*(1-exp(TauKf_ON*(t(phase_idx)-prev_end));