Why the Kf_LMax value is increased beyond its limits. Is their is any logical error Kindly help me out .

Why the Kf_LMax value is increased beyond its limits. Is their is any logical error Kindly help me out .

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 10, 5, 10, 1];  % RLC values
TauKf_ON = -0.01; % TauKf_ON
TauKf_OFF = -0.01;  % TauKf_OFF  
PhaseTimes = [0, 500, 1000, 2000, 3000, 4000, 5000]; % phase times (each row defines a phase)
% Generate time vector
t = 0:0.01:5000;
% 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', 1.5);
title('Kf_LMax over Time');
xlabel('Time');
ylabel('Kf_LMax');
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)
% Initialize output variables
Kf_LMax = zeros(size(t));

% Calculate Kf_L for each time step
for i = 1:length(t)
    Kf_LMax(i) = calculate_kf(t(i), PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);
end
% Nested function for calculating Kf_L based on time
    function Kf_L = calculate_kf(t_current, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)
        % Calculate maximum Kf_L based on RLC values using Element wise
        % division
        Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));
        
        % Default Kf_L to the maximum value of the first phase
        Kf_L = Kf_LMax_values(1);
        % Number of phases
        num_phases = numel(RLC);
        % Iterate through each phase to determine the current phase based on time
        for j = 1:num_phases
            T_start = PhaseTimes(j); % Start time of the current phase
            if j < num_phases
                T_end = PhaseTimes(j + 1); % End time of the current phase
            else
                T_end = inf; % Handle last phase separately
            end
            % Check if the current time t_current is within the current phase
            if t_current >= T_start && t_current < T_end
                if j == 1
                    % For the first phase, use the maximum value directly
                    Kf_L = Kf_LMax_values(j);
                else
                    % Time at the end of the previous phase
                    % prev_end = PhaseTimes(j - 1);
                    if j < num_phases
                        % Check RLC conditions and compute Kf_L
                        if  RLC(j - 1) < RLC(j) && RLC(j) > RLC(j + 1)
                            % Peak condition
                            Kf_endA = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) *  exp(TauKf_ON * (t_current - T_start));
                            Kf_L = Kf_LMax_values(j) + (Kf_endA - Kf_LMax_values(j - 1)) *exp(TauKf_OFF * ( t_current - T_start));
                        elseif  RLC(j - 1) < RLC(j) && RLC(j) < RLC(j + 1)
                            % Increasing RLC condition
                            Kf_endB = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t_current - T_start));
                            Kf_L = Kf_LMax_values(j) + (Kf_endB - Kf_LMax_values(j - 1)) *exp(TauKf_OFF * (t_current - T_start));
                        elseif  RLC(j - 1) > RLC(j) && RLC(j) < RLC(j + 1)
                            % Decreasing RLC condition
                            Kf_endC = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t_current - T_start));
                            Kf_L = Kf_LMax_values(j) + (Kf_endC - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t_current - T_start ));                           
                        elseif  RLC(j - 1) > RLC(j) && RLC(j) >= RLC(j + 1)
                            % Flat or decreasing RLC condition
                            Kf_endD = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t_current - T_start));
                            Kf_L = Kf_LMax_values(j) + (Kf_endD - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t_current - T_start ));        
                        end
                    else
                        % % % Handling for the last phase
                        if RLC(j - 1) < RLC(j)
                            Kf_end1 = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t_current - T_start) Kf_L = Kf_LMax_values(j) + (Kf_end1 - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t_current - T_start));
                        elseif RLC(j - 1) > RLC(j)
                            Kf_end2 = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t_current - T_start));
                            Kf_L = Kf_LMax_values(j) + (Kf_end2 - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t_current - T_start));  
                        end
                    end
                end
                return; % Exit the function after finding the correct phase
            end
        end
    end
end

Expected responses:

Answer 1

[Sam Chak]

Hi @Ehtisham

I am unfamiliar with your equation, but it appears that you did not "mathematically" impose the desired upper limit (100) on the 'Kf_L' variable. I selected a specific time at (where the peak approximately occurs) to compute the instantaneous value of 'Kf_L'​, which depends on the incremental computations of 'Kf_endA' and 'Kf_L'. The returned value is 102.6787, which clearly exceeds the desired upper limit.

In my opinion, a time-varying parameter like 'Kf_L'​, if governed by an exponentially decaying dynamic law, should be best formulated as a differential equation, and then solved either analytically or numerically.

% Define parameters
Kf_Max      = 100;  % maximum forward rate
RLC         = [0.1, 0.5, 10, 5, 10, 1];  % RLC values
TauKf_ON    = -0.01;    % TauKf_ON
TauKf_OFF   = -0.01;    % TauKf_OFF  
PhaseTimes  = [0, 500, 1000, 2000, 3000, 4000, 5000]; % phase times (each row defines a phase)

% Generate time vector
% t           = 0:0.01:5000;
t           = 1125;

% Call the function to compute receptor concentrations and Kf
[Kf_LMax]   = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)

Kf_LMax_values = 1x6

9.0909 33.3333 90.9091 83.3333 90.9091 50.0000

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

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

Kf_endA = 74.4134

Kf_L = 102.6787

Kf_LMax = 102.6787

%% TEST
% Kf_endA = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1))*exp(TauKf_ON*(t_current - T_start))
  Kf_endA = 90.9091           - (90.9091           - 33.3333              )*exp(-0.01   *(1125      - 1000))

Kf_endA = 74.4134

% Kf_L    = Kf_LMax_values(j) + (Kf_endA - Kf_LMax_values(j - 1))*exp(TauKf_OFF*(t_current - T_start));         
  Kf_L    = 90.9091           + (74.4134 - 33.3333              )*exp(-0.01    *(1125      - 1000))

Kf_L = 102.6787

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)
% Initialize output variables
Kf_LMax = zeros(size(t));

% Calculate Kf_L for each time step
for i = 1:length(t)
    Kf_LMax(i) = calculate_kf(t(i), PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);
end
% Nested function for calculating Kf_L based on time
    function Kf_L = calculate_kf(t_current, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)
        % Calculate maximum Kf_L based on RLC values using Element wise
        % division
        Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC))
        
        % Default Kf_L to the maximum value of the first phase
        Kf_L = Kf_LMax_values(1);
        
        % Number of phases
        num_phases = numel(RLC);
        % Iterate through each phase to determine the current phase based on time
        for j = 1:num_phases
            T_start = PhaseTimes(j); % Start time of the current phase
            if j < num_phases
                T_end = PhaseTimes(j + 1); % End time of the current phase
            else
                T_end = inf; % Handle last phase separately
            end
            % Check if the current time t_current is within the current phase
            if t_current >= T_start && t_current < T_end
                if j == 1
                    % For the first phase, use the maximum value directly
                    Kf_L = Kf_LMax_values(j);
                else
                    % Time at the end of the previous phase
                    % prev_end = PhaseTimes(j - 1);
                    if j < num_phases
                        % Check RLC conditions and compute Kf_L
                        if  RLC(j - 1) < RLC(j) && RLC(j) > RLC(j + 1)
                            % Peak condition
                            Kf_endA = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1))*exp(TauKf_ON * (t_current - T_start))
                            Kf_L    = Kf_LMax_values(j) + (Kf_endA - Kf_LMax_values(j - 1))*exp(TauKf_OFF*(t_current - T_start))
                        elseif  RLC(j - 1) < RLC(j) && RLC(j) < RLC(j + 1)
                            % Increasing RLC condition
                            Kf_endB = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1))*exp(TauKf_ON * (t_current - T_start));
                            Kf_L    = Kf_LMax_values(j) + (Kf_endB - Kf_LMax_values(j - 1))*exp(TauKf_OFF*(t_current - T_start));
                        elseif  RLC(j - 1) > RLC(j) && RLC(j) < RLC(j + 1)
                            % Decreasing RLC condition
                            Kf_endC = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1))*exp(TauKf_ON * (t_current - T_start));
                            Kf_L    = Kf_LMax_values(j) + (Kf_endC - Kf_LMax_values(j - 1))*exp(TauKf_OFF*(t_current - T_start ));                           
                        elseif  RLC(j - 1) > RLC(j) && RLC(j) >= RLC(j + 1)
                            % Flat or decreasing RLC condition
                            Kf_endD = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1))*exp(TauKf_ON * (t_current - T_start));
                            Kf_L    = Kf_LMax_values(j) + (Kf_endD - Kf_LMax_values(j - 1))*exp(TauKf_OFF*(t_current - T_start ));        
                        end
                    else
                        % % % Handling for the last phase
                        if RLC(j - 1) < RLC(j)
                            Kf_end1 = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t_current - T_start)); 
                            Kf_L    = Kf_LMax_values(j) + (Kf_end1 - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t_current - T_start));
                        elseif RLC(j - 1) > RLC(j)
                            Kf_end2 = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t_current - T_start));
                            Kf_L    = Kf_LMax_values(j) + (Kf_end2 - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t_current - T_start));  
                        end
                    end
                end
                return; % Exit the function after finding the correct phase
            end
        end
    end
end

[Ehtisham]

@Sam Chak Thanks for your replying .I do not want to impose the Kf_L = 100 i just want i should gives us the calculation till the Kf maximum value define in the initial condition

Answer 2

[Sam Chak]

Hi @Ehtisham

The previous answer attempts to provide a logical explanation for why the 'Kf_L'​ variable exceeds the desired upper limit (100). Are you, by any chance, looking to impose a hard limit on the 'Kf_L' variable as plotted below?

% Define parameters
Kf_Max      = 100;      % maximum forward rate
RLC         = [0.1, 0.5, 10, 5, 10, 1];  % RLC values
TauKf_ON    = -0.01;    % TauKf_ON
TauKf_OFF   = -0.01;    % TauKf_OFF  
PhaseTimes  = [0, 500, 1000, 2000, 3000, 4000, 5000]; % phase times (each row defines a phase)

% Generate time vector
t           = 0:0.01:5000;

% 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', 1.5); ylim([0 110])
title('Kf_LMax over Time');
xlabel('Time');
ylabel('Kf_LMax');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Kf_LMax] = Kf_Cal(Kf_Max, RLC, TauKf_ON, TauKf_OFF, PhaseTimes, t)
% Initialize output variables
Kf_LMax = zeros(size(t));

% Calculate Kf_L for each time step
for i = 1:length(t)
    Kf_LMax(i) = calculate_kf(t(i), PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);
end
% Nested function for calculating Kf_L based on time
    function Kf_L = calculate_kf(t_current, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)
        % Calculate maximum Kf_L based on RLC values using Element wise
        % division
        Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));
        
        % Default Kf_L to the maximum value of the first phase
        Kf_L = Kf_LMax_values(1);
        
        % Number of phases
        num_phases = numel(RLC);
        % Iterate through each phase to determine the current phase based on time
        for j = 1:num_phases
            T_start = PhaseTimes(j); % Start time of the current phase
            if j < num_phases
                T_end = PhaseTimes(j + 1); % End time of the current phase
            else
                T_end = inf; % Handle last phase separately
            end
            % Check if the current time t_current is within the current phase
            if t_current >= T_start && t_current < T_end
                if j == 1
                    % For the first phase, use the maximum value directly
                    Kf_L = Kf_LMax_values(j);
                else
                    % Time at the end of the previous phase
                    % prev_end = PhaseTimes(j - 1);
                    if j < num_phases
                        % Check RLC conditions and compute Kf_L
                        if  RLC(j - 1) < RLC(j) && RLC(j) > RLC(j + 1)
                            % Peak condition
                            Kf_endA = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1))*exp(TauKf_ON * (t_current - T_start));
                            Kf_L    = Kf_LMax_values(j) + (Kf_endA - Kf_LMax_values(j - 1))*exp(TauKf_OFF*(t_current - T_start));
                        elseif  RLC(j - 1) < RLC(j) && RLC(j) < RLC(j + 1)
                            % Increasing RLC condition
                            Kf_endB = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1))*exp(TauKf_ON * (t_current - T_start));
                            Kf_L    = Kf_LMax_values(j) + (Kf_endB - Kf_LMax_values(j - 1))*exp(TauKf_OFF*(t_current - T_start));
                        elseif  RLC(j - 1) > RLC(j) && RLC(j) < RLC(j + 1)
                            % Decreasing RLC condition
                            Kf_endC = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1))*exp(TauKf_ON * (t_current - T_start));
                            Kf_L    = Kf_LMax_values(j) + (Kf_endC - Kf_LMax_values(j - 1))*exp(TauKf_OFF*(t_current - T_start ));                           
                        elseif  RLC(j - 1) > RLC(j) && RLC(j) >= RLC(j + 1)
                            % Flat or decreasing RLC condition
                            Kf_endD = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1))*exp(TauKf_ON * (t_current - T_start));
                            Kf_L    = Kf_LMax_values(j) + (Kf_endD - Kf_LMax_values(j - 1))*exp(TauKf_OFF*(t_current - T_start ));        
                        end
                    else
                        % % % Handling for the last phase
                        if RLC(j - 1) < RLC(j)
                            Kf_end1 = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t_current - T_start)); 
                            Kf_L    = Kf_LMax_values(j) + (Kf_end1 - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t_current - T_start));
                        elseif RLC(j - 1) > RLC(j)
                            Kf_end2 = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t_current - T_start));
                            Kf_L    = Kf_LMax_values(j) + (Kf_end2 - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t_current - T_start));  
                        end
                    end
                end
                Kf_L    = min(Kf_L, Kf_Max);    % Sam: imposing a hard limit Kf_Max
                return; % Exit the function after finding the correct phase
            end
        end
    end
end

[Ehtisham]

@Sam Chak I just need help that is their is any logical error in the code beacuse if even we vary the Kf_Max but the output is crossing the limited values. Kindly guide me Thanks

[Sam Chak]

@Ehtisham, Could you provide a sketch of the expected profile of 'Kf_L'?

[Torsten]

@Ehtisham

How should we be able to decide whether your formula to compute Kf_L is correct or not ?

The usual exponential between an initial state u0 and a final state uinf is

u(t) = uinf + (u0-uinf)*exp(-c*(t-t0)) 
c>0 a numerical factor

and u(t) will always run between unif and u0 for t>=t0.

[Ehtisham]

@Sam Chak

Now i used this xponentially decaying dynamic law, formulated as a differential equation, but i am not getting the dynamic behaviour when i change the values of RLC

[Sam Chak]

Hi @Ehtisham

I’ve put your code here so that other users can check it out.

Have you been able to find an explanation for the responses from a purely mathematical perspective?

By the way, you haven’t provided the requested sketch of the expected responses.

%% Code in SimfileNPhaseRun.m
% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 10, 1, 0.2, 0.4];  % RLC values
TauKf_ON = -0.1; % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF
Kb_Max = 50;  % maximum backward rate  
PhaseTimes = [0, 500, 1000, 1500, 2000, 3000, 4000];% phase times (each row defines a phase)
timespan = [0,  4000];  % 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)

[t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase (Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...
    Kb_Max,  PhaseTimes, timespan, Non_Active_Receptor_concentration, Active_Receptor_concentration);

%% Code in SimfileNPhase.m 
function [t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase (Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...
    Kb_Max,   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(Kb_Max);

    % 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Nested function for calculating Kf_L based on time
    function Kf_L = calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)
        % Calculate maximum Kf_L based on RLC values using element-wise division
        Kf_LMax_values = Kf_Max * (RLC ./ (1 + RLC));
        
        % Default Kf_L to the maximum value of the first phase
        Kf_L = Kf_LMax_values(1);

        % Number of phases
        num_phases = numel(RLC);

        % Iterate through each phase to determine the current phase based on time
        for j = 1:num_phases
            T_start = PhaseTimes(j); % Start time of the current phase
            if j < num_phases
                T_end = PhaseTimes(j + 1); % End time of the current phase
            else
                T_end = inf; % Handle last phase separately
            end

            % Check if the current time t_current is within the current phase
            if t >= T_start && t < T_end
                if j == 1
                    % For the first phase, use the maximum value directly
                    Kf_L = Kf_LMax_values(j);
                else
                    % Time at the end of the previous phase
                    prev_end = PhaseTimes(j - 1);
                    if j < num_phases
                        % Check RLC conditions and compute Kf_L
                        if RLC(j - 1) < RLC(j) && RLC(j) > RLC(j + 1)
                            % Peak condition
                            Kf_endA = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t - prev_end));
                            Kf_L = Kf_LMax_values(j) + (Kf_endA - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t - prev_end));
                        elseif RLC(j - 1) < RLC(j) && RLC(j) < RLC(j + 1)
                            % Increasing RLC condition
                            Kf_endB = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t - prev_end));
                            Kf_L = Kf_LMax_values(j) + (Kf_endB - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t - prev_end));
                        elseif RLC(j - 1) > RLC(j) && RLC(j) < RLC(j + 1)
                            % Decreasing RLC condition
                            Kf_endC = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t - prev_end));
                            Kf_L = Kf_LMax_values(j) + (Kf_endC - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t - prev_end ));                           
                        elseif RLC(j - 1) > RLC(j) && RLC(j) >= RLC(j + 1)
                            % Flat or decreasing RLC condition
                            Kf_endD = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t - prev_end));
                            Kf_L = Kf_LMax_values(j) + (Kf_endD - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t - prev_end));        
                        end
                     else
                        % % % Handling for the last phase
                        if RLC(j - 1) < RLC(j)
                            Kf_end1 = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t - prev_end));
                             Kf_L = Kf_LMax_values(j) + (Kf_end1 - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t - prev_end));
                        elseif RLC(j - 1) > RLC(j)
                            Kf_end2 = Kf_LMax_values(j) - (Kf_LMax_values(j) - Kf_LMax_values(j - 1)) * exp(TauKf_ON * (t - prev_end));
                            Kf_L = Kf_LMax_values(j) + (Kf_end2 - Kf_LMax_values(j - 1)) * exp(TauKf_OFF * (t - prev_end));  
                        end
                    end
                end

                return; % Exit the function after finding the correct phase
            
            end
           
        end
    end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kb = calculate_kb(Kb_Max)
  
Kb = Kb_Max;  % Default to the maximum value
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]

here the skecth of active and non active receptor and reaction rates @Sam Chak

[Sam Chak]

Hi @Ehtisham

I’m not familiar with the concentration system. However, since this is a super simple coupled first-order system, you can probably run the simulation phase by phase—integrate, stop, and then integrate again until the next phase.

Kf_L    = 1.5e-3;
Kb      = Kf_L*9;
tspan   = [0, 500];
y0      = [100; 0];
[t, y]  = ode45(@(t, y) ode_LR(t, y, Kf_L, Kb), tspan, y0);

tl  = tiledlayout(2, 1, 'TileSpacing', 'Compact');
nexttile
plot(t, y(:,1)), grid on
title('Non-Active Receptor concentration')
nexttile
plot(t, y(:,2)), grid on
title('Active Receptor concentration')
xlabel(tl, 'Time')

function dydt = ode_LR(t, y, Kf_L, Kb)
    % Non_Active_Receptor_concentration   = y(1);
    % Active_Receptor_concentration       = y(2);
    dy1     = - Kf_L*y(1) + Kb*y(2);
    dy2     =   Kf_L*y(1) - Kb*y(2);
    dydt    = [dy1; 
               dy2];
end

Answer 3

[Sam Chak]

Hi @Ehtisham

If the gains and stay constant or vary very slowly during each phase, then you can consider applying the following analytical solutions by making the appropriate substitutions for the gains and initial conditions at each phase.

syms x(t) y(t) K1 K2 C1 C2 

% x  is the concentration of Non-Active Receptor
% y  is the concentration of Active Receptor
% K1 is Kf_L
% K2 is Kb
% C1 is initial condition of x
% C2 is initial condition of y

eqns = [diff(x,t) == - K1*x + K2*y,
        diff(y,t) ==   K1*x - K2*y];

cond = [x(0)==C1, 
        y(0)==C2];

[xSol(t), ySol(t)] = dsolve(eqns, cond)

xSol(t) = 

ySol(t) =