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

Sign in to answer this question.

Follow Question

Answer 1

Sam Chak on 10 Sep 2024

0 votes

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

1 Comment
Show -1 older comments Hide -1 older comments

Ehtisham on 10 Sep 2024

@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 on 10 Sep 2024

0 votes

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

Kf_L = min(Kf_L, Kf_Max); % Sam: imposing a hard limit Kf_Max

return; % Exit the function after finding the correct phase

end

7 Comments
Show 5 older comments Hide 5 older comments

Sam Chak on 11 Sep 2024

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

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

% 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

return; % Exit the function after finding the correct phase

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

Ehtisham on 11 Sep 2024

Edited: Sam Chak on 11 Sep 2024

IMG_20240911_134411_831.jpg

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

Sam Chak on 11 Sep 2024

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 on 11 Sep 2024

0 votes

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) =

0 Comments
Show -2 older comments Hide -2 older comments