When i call and run this code it just save the Phase1 results in CSvV file not other results .
85 ビュー (過去 30 日間)
古いコメントを表示
When i call and run this code it just save the Phase1 results in CSvV file not other results .
2 件のコメント
Jaimin
2024 年 8 月 12 日 10:26
Could you please provide detailed information about the expected results in CSV format? Thank you.
採用された回答
Neelanshu
2024 年 8 月 12 日 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 件のコメント
Neelanshu
2024 年 8 月 13 日 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!
その他の回答 (1 件)
Torsten
2024 年 8 月 12 日 10:11
編集済み: Torsten
2024 年 8 月 12 日 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 件のコメント
Torsten
2024 年 8 月 12 日 12:45
編集済み: Torsten
2024 年 8 月 12 日 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
参考
カテゴリ
Help Center および File Exchange で Image Processing and Computer Vision についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!