When i call and run this code it just save the Phase1 results in CSvV file not other results .
85 views (last 30 days)
Show older comments
Ehtisham
on 12 Aug 2024 at 9:51
Commented: Ehtisham
on 14 Aug 2024 at 8:57
When i call and run this code it just save the Phase1 results in CSvV file not other results .
2 Comments
Jaimin
on 12 Aug 2024 at 10:26
Could you please provide detailed information about the expected results in CSV format? Thank you.
Accepted Answer
Neelanshu
on 12 Aug 2024 at 10:28
Hi Ehtisham,
I understand that you are unable to save CSV results other than 'Phase1.csv'.
The issue occurs because the PhaseResults variable is a 1x1 struct. This happened as the array used in 'arrayfun', as shown in the following code from TestCode1.m, contains only 1 element :
% Calculate phase results
PhaseResults = arrayfun(@(i) ...
calculate_phase(t, Active_Receptor_concentration, PhaseTimes(i, :), RLC(i)) ...
, 1:size(PhaseTimes, 1), 'UniformOutput', false);
The size of PhaseTimes, a row vector, is [1 8] which leads to the array utilized in the 'arrayfun' to have a single element. To resolve the issue, kindly ensure that the array of the correct dimension is created.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1751054/image.png)
You may refer to the following documentation to learn more about 'arrayfun' function:
Hope this helps!
4 Comments
Neelanshu
on 13 Aug 2024 at 9:40
Sure I will elaborate the cause of error. The issues occurs because the 'calculate_phase' function is not compatible with row vectors like 'PhaseTimes'.
I have updated the code at the following sections :
PhaseResults = arrayfun(@(i) ...
calculate_phase(t, Active_Receptor_concentration, PhaseTimes(:,i:i+1), RLC(i)) ...
, 1:(size(PhaseTimes, 2)-1), 'UniformOutput', false);
This sections creates a 1x7 array and calculate_phase works on the elements of this array.
function result = calculate_phase(t, Active_Receptor_concentration, PhaseTimes, RLC)
T_start = PhaseTimes(:,1);
T_end = PhaseTimes(:,2);
In this section, calculate_phase is modified to work with row vectors.
This results in the formation of 7 CSV files in the current path as shown in the image :
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1751599/image.png)
Here is the entire code script:
% Define parameters
Kf_Max = 100; % maximum forward rate
RLC = [0.01, 1, 5, 1, 7, 10, 0.01]; % RLC values
TauKf_ON = -0.01; % TauKf_ON
TauKf_OFF = -0.01; % TauKf_OFF
Kb_Max = 80; % maximum backward rate
Kb_Min = 10; % minimum backward rate
TauKb_ON = -0.01; % TauKb_ON
TauKb_OFF = -0.01; % TauKb_OFF
PhaseTimes = [0, 10, 500, 1000, 2000, 3000, 4000,5000];% phase times (each row defines a phase)
timespan = [0, 5000]; % timespan for the simulation
% Example initial conditions for non-active and active receptor concentrations
Non_Active_Receptor_concentration = 100; % Initial non-active receptor concentration (as a vector)
Active_Receptor_concentration =0; % Initial active receptor concentration (as a vector)
% Run the function
[t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = TestCode1(Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...
Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, PhaseTimes, timespan, Non_Active_Receptor_concentration, Active_Receptor_concentration);
function [t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = TestCode1(Kf_Max, RLC, TauKf_ON, ~, ...
Kb_Max, Kb_Min, TauKb_ON, ~, PhaseTimes, timespan, Non_Active_Receptor_concentration, Active_Receptor_concentration)
% Define functions for forward and backward reaction rates
Kf_L = @(t) calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON);
Kb = @(t) calculate_kb(t, PhaseTimes, Kb_Max, Kb_Min, TauKb_ON, RLC);
% Ensure that Non_Active_Receptor_concentration and Active_Receptor_concentration are column vectors
Non_Active_Receptor_concentration = Non_Active_Receptor_concentration(:);
Active_Receptor_concentration = Active_Receptor_concentration(:);
% Set initial conditions
initial_conditions = [Non_Active_Receptor_concentration(1); Active_Receptor_concentration(1)];
% Set ODE options for step size
options = odeset('MaxStep', 0.05, 'RelTol', 1e-6, 'AbsTol', 1e-8);
% Solve the ODE system
[t, y] = ode45(@(t, y) ode_LR(t, y, Kf_L, Kb), timespan, initial_conditions, options);
% Extract the concentrations
Non_Active_Receptor_concentration = y(:, 1);
Active_Receptor_concentration = y(:, 2);
% Calculate forward and backward reaction rates over time
Kf_values = arrayfun(Kf_L, t);
Kb_values = arrayfun(Kb, t);
% Plot active and non-active receptor concentrations
figure;
plot(t, Non_Active_Receptor_concentration, 'r', 'DisplayName', 'Non-Active Receptor Concentration');
hold on;
plot(t, Active_Receptor_concentration, 'b', 'DisplayName', 'Active Receptor Concentration');
legend;
xlabel('Time');
ylabel('Concentration');
title('Receptor Concentrations');
hold off;
% Plot forward and backward reaction rates
figure;
plot(t, Kf_values, 'k', 'DisplayName', 'Forward Reaction Rate');
hold on;
plot(t, Kb_values, 'c', 'DisplayName', 'Backward Reaction Rate');
legend;
xlabel('Time');
ylabel('Reaction Rate');
title('Reaction Rates');
hold off;
% Calculate phase results
PhaseResults = arrayfun(@(i) ...
calculate_phase(t, Active_Receptor_concentration, PhaseTimes(:,i:i+1), RLC(i)) ...
, 1:(size(PhaseTimes, 2)-1), 'UniformOutput', false);
PhaseResults = vertcat(PhaseResults{:}); % Concatenate results
% Calculate maximum forward reaction rate
Kf_LMax = Kf_Max * (RLC ./ (RLC + 1));
% Write Phase results to CSV
for i = 1:size(PhaseResults, 1)
PhaseTable = struct2table(PhaseResults(i))
writetable(PhaseTable, ['Phase', num2str(i), '.csv']);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kf_L = calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON)
Kf_LMax = Kf_Max * (RLC ./ (RLC + 1));
Kf_L = Kf_LMax(1); % Default to the first phase value
num_phases = numel(RLC);
for j = 1:num_phases
T_start = PhaseTimes(j);
T_end = PhaseTimes(j + 1);
if t >= T_start && t < T_end
if j == 1
Kf_L = Kf_LMax(j);
else
prev_end = PhaseTimes(j);
if j < num_phases
if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
Kf_L = Kf_LMax(j);
elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
end
end
end
return;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kb = calculate_kb(t, PhaseTimes, Kb_Max, Kb_Min, TauKb_ON, RLC)
Kb = Kb_Min; % Default to the minimum value
if all(RLC == RLC(1))
Kb = Kb_Min; % Set Kb to Kb_Min if all RLC values are equal
return;
end
for j = 1:numel(RLC)
T_start = PhaseTimes(j);
T_end = PhaseTimes(j + 1);
if t >= T_start && t < T_end
if j == 1
Kb = Kb_Min;
else
% prev_end = PhaseTimes(j);
if j < numel(RLC)
if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
end
else
if RLC(j-1) < RLC(j)
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
elseif RLC(j-1) > RLC(j)
Kb = Kb_Max - (Kb_Max - Kb_Min ) * exp(TauKb_ON * (T_end - T_start));
end
end
end
return;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dydt = ode_LR(t, y, Kf_L, Kb)
Non_Active_Receptor_concentration = y(1);
Active_Receptor_concentration = y(2);
dNonActiveReceptor_dt = -Kf_L(t) * Non_Active_Receptor_concentration + Kb(t) * Active_Receptor_concentration;
dActiveReceptor_dt = Kf_L(t) * Non_Active_Receptor_concentration - Kb(t) * Active_Receptor_concentration;
dydt = [dNonActiveReceptor_dt; dActiveReceptor_dt];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result = calculate_phase(t, Active_Receptor_concentration, PhaseTimes, RLC)
T_start = PhaseTimes(:,1);
T_end = PhaseTimes(:,2);
Phase_indices = (t >= T_start & t <= T_end);
Phase_concentration = Active_Receptor_concentration(Phase_indices);
Phase_time = t(Phase_indices);
% Calculate peak and steady-state values
[Rpeak, Tpeak, peak_type] = findpeak(Phase_time, Phase_concentration);
Rss = mean(Active_Receptor_concentration(t >= (T_end - 10) & t <= T_end));
% Calculate the T50 (time to reach half of peak value)
half_peak_value = (Rss + Rpeak) / 2;
[~, idx_50_percent] = min(abs(Phase_concentration - half_peak_value));
T50 = Phase_time(idx_50_percent) - Tpeak;
% Calculate other metrics
ratio_Rpeak_Rss = Rpeak / Rss;
Delta = Rpeak - Rss;
L = RLC;
% Package results in a struct
result.Rpeak = Rpeak;
result.Rss = Rss;
result.Tpeak = Tpeak;
result.T50 = T50;
result.ratio_Rpeak_Rss = ratio_Rpeak_Rss;
result.Delta = Delta;
result.L = L;
result.peak_type = peak_type;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Rpeak, Tpeak, peak_type] = findpeak(time, concentration)
% Compute the derivative
dt = diff(time);
dy = diff(concentration);
derivative = dy ./ dt;
% Find zero-crossings of the derivative
zero_crossings = find(diff(sign(derivative)) ~= 0);
% Initialize output variables
Rpeak = NaN;
Tpeak = NaN;
peak_type = 'None';
% Check if there are zero crossings
if ~isempty(zero_crossings)
% Identify peaks and troughs by examining the sign changes
max_indices = zero_crossings(derivative(zero_crossings) > 0 & derivative(zero_crossings + 1) < 0);
min_indices = zero_crossings(derivative(zero_crossings) < 0 & derivative(zero_crossings + 1) > 0);
% Check if there are maxima or minima
if ~isempty(max_indices) || ~isempty(min_indices)
if ~isempty(max_indices) && ~isempty(min_indices)
% Find the highest maximum
[Rpeak_max, maxIdx] = max(concentration(max_indices));
% Find the lowest minimum
[Rpeak_min, minIdx] = min(concentration(min_indices));
% Determine whether the highest maximum is greater or the lowest minimum is smaller
if Rpeak_max >= abs(Rpeak_min)
Rpeak = Rpeak_max;
Tpeak = time(max_indices(maxIdx));
peak_type = 'High';
else
Rpeak = Rpeak_min;
Tpeak = time(min_indices(minIdx));
peak_type = 'Low';
end
% If only maxima exist
elseif ~isempty(max_indices)
[Rpeak, maxIdx] = max(concentration(max_indices));
Tpeak = time(max_indices(maxIdx));
peak_type = 'High';
% If only minima exist
elseif ~isempty(min_indices)
[Rpeak, minIdx] = min(concentration(min_indices));
Tpeak = time(min_indices(minIdx));
peak_type = 'Low';
end
end
end
end
Hope this helps!
More Answers (1)
Torsten
on 12 Aug 2024 at 10:11
Edited: Torsten
on 12 Aug 2024 at 10:12
Use
PhaseTable = struct2table(PhaseResults);
writetable(PhaseTable,'Phase.csv');
instead of
% Write Phase results to CSV
for i = 1:size(PhaseResults, 1)
PhaseTable = struct2table(PhaseResults(i));
writetable(PhaseTable, ['Phase', num2str(i), '.csv']);
end
You only got one struct with name "PhaseResults".
And if you had several structs PhaseResults, you couldn't access them by "PhaseResults(i)", but by "PhaseResults{i}".
3 Comments
Torsten
on 12 Aug 2024 at 12:45
Edited: Torsten
on 12 Aug 2024 at 12:58
% Define parameters
Kf_Max = 100; % maximum forward rate
RLC = [0.01, 1, 5, 1, 7, 10, 0.01]; % RLC values
TauKf_ON = -0.01; % TauKf_ON
TauKf_OFF = -0.01; % TauKf_OFF
Kb_Max = 80; % maximum backward rate
Kb_Min = 10; % minimum backward rate
TauKb_ON = -0.01; % TauKb_ON
TauKb_OFF = -0.01; % TauKb_OFF
PhaseTimes = [0, 10, 500, 1000, 2000, 3000, 4000,5000];% phase times (each row defines a phase)
timespan = [0, 5000]; % timespan for the simulation
% Example initial conditions for non-active and active receptor concentrations
Non_Active_Receptor_concentration = 100; % Initial non-active receptor concentration (as a vector)
Active_Receptor_concentration =0; % Initial active receptor concentration (as a vector)
% Run the function
[t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = TestCode1(Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...
Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, PhaseTimes, timespan, Non_Active_Receptor_concentration, Active_Receptor_concentration);
function [t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = TestCode1(Kf_Max, RLC, TauKf_ON, ~, ...
Kb_Max, Kb_Min, TauKb_ON, ~, PhaseTimes, timespan, Non_Active_Receptor_concentration, Active_Receptor_concentration)
% Define functions for forward and backward reaction rates
Kf_L = @(t) calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON);
Kb = @(t) calculate_kb(t, PhaseTimes, Kb_Max, Kb_Min, TauKb_ON, RLC);
% Ensure that Non_Active_Receptor_concentration and Active_Receptor_concentration are column vectors
Non_Active_Receptor_concentration = Non_Active_Receptor_concentration(:);
Active_Receptor_concentration = Active_Receptor_concentration(:);
% Set initial conditions
initial_conditions = [Non_Active_Receptor_concentration(1); Active_Receptor_concentration(1)];
% Set ODE options for step size
options = odeset('MaxStep', 0.05, 'RelTol', 1e-6, 'AbsTol', 1e-8);
% Solve the ODE system
[t, y] = ode45(@(t, y) ode_LR(t, y, Kf_L, Kb), timespan, initial_conditions, options);
% Extract the concentrations
Non_Active_Receptor_concentration = y(:, 1);
Active_Receptor_concentration = y(:, 2);
% Calculate forward and backward reaction rates over time
Kf_values = arrayfun(Kf_L, t);
Kb_values = arrayfun(Kb, t);
% Plot active and non-active receptor concentrations
figure;
plot(t, Non_Active_Receptor_concentration, 'r', 'DisplayName', 'Non-Active Receptor Concentration');
hold on;
plot(t, Active_Receptor_concentration, 'b', 'DisplayName', 'Active Receptor Concentration');
legend;
xlabel('Time');
ylabel('Concentration');
title('Receptor Concentrations');
hold off;
% Plot forward and backward reaction rates
figure;
plot(t, Kf_values, 'k', 'DisplayName', 'Forward Reaction Rate');
hold on;
plot(t, Kb_values, 'c', 'DisplayName', 'Backward Reaction Rate');
legend;
xlabel('Time');
ylabel('Reaction Rate');
title('Reaction Rates');
hold off;
% Calculate maximum forward reaction rate
Kf_LMax = Kf_Max * (RLC ./ (RLC + 1));
for j = 1:size(PhaseTimes,1)
% Calculate phase results
for i = 1:size(PhaseTimes,2)-1
PhaseResults{i} = calculate_phase(t, Active_Receptor_concentration, PhaseTimes(j,i:i+1), RLC(i));
end
PhaseResults = vertcat(PhaseResults{:}); % Concatenate results
% Write Phase results to CSV
PhaseTable = struct2table(PhaseResults)
writetable(PhaseTable, ['Phase', num2str(j), '.csv']);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kf_L = calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON)
Kf_LMax = Kf_Max * (RLC ./ (RLC + 1));
Kf_L = Kf_LMax(1); % Default to the first phase value
num_phases = numel(RLC);
for j = 1:num_phases
T_start = PhaseTimes(j);
T_end = PhaseTimes(j + 1);
if t >= T_start && t < T_end
if j == 1
Kf_L = Kf_LMax(j);
else
prev_end = PhaseTimes(j);
if j < num_phases
if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
Kf_L = Kf_LMax(j);
elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j-1)) * exp(TauKf_ON * (t - prev_end));
end
end
end
return;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kb = calculate_kb(t, PhaseTimes, Kb_Max, Kb_Min, TauKb_ON, RLC)
Kb = Kb_Min; % Default to the minimum value
if all(RLC == RLC(1))
Kb = Kb_Min; % Set Kb to Kb_Min if all RLC values are equal
return;
end
for j = 1:numel(RLC)
T_start = PhaseTimes(j);
T_end = PhaseTimes(j + 1);
if t >= T_start && t < T_end
if j == 1
Kb = Kb_Min;
else
% prev_end = PhaseTimes(j);
if j < numel(RLC)
if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
end
else
if RLC(j-1) < RLC(j)
Kb = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (T_end - T_start));
elseif RLC(j-1) > RLC(j)
Kb = Kb_Max - (Kb_Max - Kb_Min ) * exp(TauKb_ON * (T_end - T_start));
end
end
end
return;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dydt = ode_LR(t, y, Kf_L, Kb)
Non_Active_Receptor_concentration = y(1);
Active_Receptor_concentration = y(2);
dNonActiveReceptor_dt = -Kf_L(t) * Non_Active_Receptor_concentration + Kb(t) * Active_Receptor_concentration;
dActiveReceptor_dt = Kf_L(t) * Non_Active_Receptor_concentration - Kb(t) * Active_Receptor_concentration;
dydt = [dNonActiveReceptor_dt; dActiveReceptor_dt];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result = calculate_phase(t, Active_Receptor_concentration, PhaseTimes, RLC)
T_start = PhaseTimes(1);
T_end = PhaseTimes(2);
Phase_indices = (t >= T_start & t <= T_end);
Phase_concentration = Active_Receptor_concentration(Phase_indices);
Phase_time = t(Phase_indices);
% Calculate peak and steady-state values
[Rpeak, Tpeak, peak_type] = findpeak(Phase_time, Phase_concentration);
Rss = mean(Active_Receptor_concentration(t >= (T_end - 10) & t <= T_end));
% Calculate the T50 (time to reach half of peak value)
half_peak_value = (Rss + Rpeak) / 2;
[~, idx_50_percent] = min(abs(Phase_concentration - half_peak_value));
T50 = Phase_time(idx_50_percent) - Tpeak;
% Calculate other metrics
ratio_Rpeak_Rss = Rpeak / Rss;
Delta = Rpeak - Rss;
L = RLC;
% Package results in a struct
result.Rpeak = Rpeak;
result.Rss = Rss;
result.Tpeak = Tpeak;
result.T50 = T50;
result.ratio_Rpeak_Rss = ratio_Rpeak_Rss;
result.Delta = Delta;
result.L = L;
result.peak_type = peak_type;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Rpeak, Tpeak, peak_type] = findpeak(time, concentration)
% Compute the derivative
dt = diff(time);
dy = diff(concentration);
derivative = dy ./ dt;
% Find zero-crossings of the derivative
zero_crossings = find(diff(sign(derivative)) ~= 0);
% Initialize output variables
Rpeak = NaN;
Tpeak = NaN;
peak_type = 'None';
% Check if there are zero crossings
if ~isempty(zero_crossings)
% Identify peaks and troughs by examining the sign changes
max_indices = zero_crossings(derivative(zero_crossings) > 0 & derivative(zero_crossings + 1) < 0);
min_indices = zero_crossings(derivative(zero_crossings) < 0 & derivative(zero_crossings + 1) > 0);
% Check if there are maxima or minima
if ~isempty(max_indices) || ~isempty(min_indices)
if ~isempty(max_indices) && ~isempty(min_indices)
% Find the highest maximum
[Rpeak_max, maxIdx] = max(concentration(max_indices));
% Find the lowest minimum
[Rpeak_min, minIdx] = min(concentration(min_indices));
% Determine whether the highest maximum is greater or the lowest minimum is smaller
if Rpeak_max >= abs(Rpeak_min)
Rpeak = Rpeak_max;
Tpeak = time(max_indices(maxIdx));
peak_type = 'High';
else
Rpeak = Rpeak_min;
Tpeak = time(min_indices(minIdx));
peak_type = 'Low';
end
% If only maxima exist
elseif ~isempty(max_indices)
[Rpeak, maxIdx] = max(concentration(max_indices));
Tpeak = time(max_indices(maxIdx));
peak_type = 'High';
% If only minima exist
elseif ~isempty(min_indices)
[Rpeak, minIdx] = min(concentration(min_indices));
Tpeak = time(min_indices(minIdx));
peak_type = 'Low';
end
end
end
end
See Also
Categories
Find more on Image Processing and Computer Vision in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!