(Kf_LMax_values(i) - Kf_LMax_values(i - 1)) --> 0 because you set the fisrt values all to Kf_LMax_values
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?
135 次查看（过去 30 天）
显示 更早的评论
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
5 个评论
Ehtisham
2024-10-8
编辑：dpb
2024-10-8
@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
Voss
2024-10-8
"(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
2024-10-10
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');
采纳的回答
Voss
2024-10-8
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
27 个评论
Ehtisham
2024-10-11
@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
2024-10-11
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
2024-10-11
@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
2024-10-14
编辑：Torsten
2024-10-14
@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
2024-10-14
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
2024-10-16
编辑：Voss
2024-10-16
@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
2024-10-16
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.
Torsten
2024-10-21，9:45
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
2024-10-21，12:28
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
2024-10-22，9:49
@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
2024-10-22，11:54
编辑：Torsten
2024-10-22，11:55
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
2024-10-22，12:05
编辑：Ehtisham
2024-10-22，12:06
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
2024-10-22，12:24
编辑：Torsten
2024-10-23，0:28
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
2024-10-24，6:42
编辑：Walter Roberson
2024-10-30，18:40
@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
2024-10-30，18:44
But I gave you the answer: Because "last_phase_idx" becomes empty in the course of your calculations.
Walter Roberson
2024-11-1，7:35
After
phase_idx = (t >= T_start) & (t < T_end);
insert
if ~any(phase_idx); continue; end
Ehtisham
2024-11-1，7:59
编辑：Walter Roberson
2024-11-7，18:35
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
2024-11-7，18:48
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.
Torsten
2024-11-11，11:43
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
2024-11-11，17:44
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
2024-11-12，8:03
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
更多回答（1 个）
dpb
2024-10-8
移动：dpb
2024-10-8
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));
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Plot Customization 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)