![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1596906/image.png)
getting error in code when running
2 ビュー (過去 30 日間)
古いコメントを表示
CheckKFmax_1()
function CheckKFmax_1()
% Parameters
timespan = 0:0.1:500; % Time vector
Receptor_concentration = 100; % Concentration of Receptor
% Initial conditions
C_LigandReceptor_0 = 0; % Initial concentration of the ligand-receptor complex
% Define different ranges of Kf1Max values
Kf1Max_values = [20, 40, 60, 80, 100];
Kb1Max_values = [20, 40, 60, 80, 100];
% Create a cell array to store tables for different Kf1Max and Kb1Max ranges
tables_cell = cell(length(Kf1Max_values), length(Kb1Max_values));
for i = 1:length(Kf1Max_values)
Kf1Max = Kf1Max_values(i);
for j = 1:length(Kb1Max_values)
Kb1Max = Kb1Max_values(j);
% Time-dependent forward reaction rate constants
t_start = 5; % Start time of receptor pulse
t_end = 500; % Duration of receptor pulse
Kf1Min = 1;
Kf1Tau_on = -1;
Kf1Tau_off = -1;
% Time-dependent backward reaction rate constants
Kb1Min = 10;
Kb1Tau_on = -0.01;
Kb1Tau_off = -0.01;
% Define anonymous functions for forward and backward reaction rates
kf_1 = @(t) calculate_kf(t, t_start, t_end, Kf1Max, Kf1Min, Kf1Tau_on, Kf1Tau_off);
kb_1 = @(t) calculate_kb(t, t_start, t_end, Kb1Max, Kb1Min, Kb1Tau_on, Kb1Tau_off);
% Solve the ODE system
[t, y] = ode45(@(t, y) ode_LR(t, y, kf_1, kb_1), timespan, [Receptor_concentration; C_LigandReceptor_0]);
% Extract the concentrations
receptor_concentration = y(:, 1);
ligand_Receptor_concentration = y(:, 2);
% Displaying the peak value and its corresponding time
[peak_value, peak_time_idx] = max(ligand_Receptor_concentration);
time_to_peak = t(peak_time_idx);
% Calculate the time between start and peak time
time_to_peak = time_to_peak - t_start;
% Find steady state and time to approach 50% of the steady state
steady_state = mean(ligand_Receptor_concentration(end-50:end)); % Considering the last few points for steady state
% Calculate T-50: Time to reach 50% of the peak value
half_peak_value = (steady_state + peak_value)/2;
% Find the indices where the ligand receptor concentration is closest to half the peak value
[~, idx_50_percent] = min(abs(ligand_Receptor_concentration(peak_time_idx:end) - half_peak_value));
% Get the time corresponding to the closest index
time_to_50_percent = t(idx_50_percent)- time_to_peak;
% Find the ratio of the R*peak/R* steady state
peak_to_steady_state_ratio = peak_value / steady_state;
% Find the ratio of the R*peak/R* steady state
Delta = peak_value - steady_state;
disp(['R*peak: ', num2str(peak_value)]);
disp(['T-peak: ', num2str(time_to_peak)]);
disp(['50% form peak to ss : ', num2str(half_peak_value)]);
disp(['R* Steady state: ', num2str(steady_state)]);
disp(['T-50 : ', num2str(time_to_50_percent)]);
disp(['R*peak/R*ss: ', num2str(peak_to_steady_state_ratio)]);
disp(['Δ (R*peak-R*ss) : ', num2str(Delta)]);
R_peak_values(i) = peak_value;
T_peak_values(i) = time_to_peak;
R_steady_state_values(i) = steady_state;
T_50_values(i) = time_to_50_percent;
R_peak_R_ss_values(i) = peak_to_steady_state_ratio;
Delta_values(i) = Delta;
% Create a table to combine all results
data_table = table(R_peak_values', T_peak_values', R_steady_state_values', T_50_values', R_peak_R_ss_values', Delta_values', ...
'VariableNames', {'R_peak', 'T_peak', 'R_steady_state', 'T_50', 'R_peak_R_ss', 'Delta'});
% Store the table in the cell array for each Kf1Max range
tables_cell{i, j} = data_table;
% Generate a unique filename for each combination of Kf1Max and Kb1Max values
csv_filename = sprintf('interaction_data_Kf1Max_%d_Kb1Max_%d.csv', Kf1Max, Kb1Max);
writetable(data_table, csv_filename);
end
end
% Combine all results into a single table
combined_table = vertcat(tables_cell{:});
% Generate a unique filename for the combined table
combined_csv_filename = 'combined_interaction_data_all.csv';
writetable(combined_table, combined_csv_filename);
% Load the combined interaction data CSV file
combined_csv_filename = 'combined_interaction_data_all.csv';
combined_table = readtable(combined_csv_filename);
% Extract Kf1Max and Kb1Max values from the combined table
Kf1Max_values = unique(combined_table.Kf1Max);
Kb1Max_values = unique(combined_table.Kb1Max);
% Initialize matrices to store Kf1Max, Kb1Max, and R*peak values
Kf1Max_matrix = zeros(length(Kb1Max_values), length(Kf1Max_values));
Kb1Max_matrix = zeros(length(Kb1Max_values), length(Kf1Max_values));
R_peak_matrix = zeros(length(Kb1Max_values), length(Kf1Max_values));
% Populate matrices with values from the combined table
for i = 1:length(Kf1Max_values)
for j = 1:length(Kb1Max_values)
idx = find(combined_table.Kf1Max == Kf1Max_values(i) & combined_table.Kb1Max == Kb1Max_values(j));
if ~isempty(idx)
Kf1Max_matrix(j, i) = Kf1Max_values(i);
Kb1Max_matrix(j, i) = Kb1Max_values(j);
R_peak_matrix(j, i) = combined_table.R_peak(idx);
end
end
end
% Create a 3D fsurf plot
figure;
fsurf(Kf1Max_matrix, Kb1Max_matrix, R_peak_matrix);
xlabel('Kf1Max');
ylabel('Kb1Max');
zlabel('R*peak');
title('3D fsurf Plot for Combined R*peak Values');
end
function kf_1 = calculate_kf(t, t_start, t_end, Kf1Max, Kf1Min, Kf1Tau_on, Kf1Tau_off)
% Initialize kf_1
kf_1 = zeros(size(t));
% Initial phase: Constant value Kf1Min
initial_phase = t < t_start;
kf_1(initial_phase) = Kf1Min;
% Increasing phase: Exponential growth from Kf1Min to Kf1Max
increasing_phase = t >= t_start & t < t_end;
kf_1(increasing_phase) = Kf1Max - (Kf1Max - Kf1Min) * exp(Kf1Tau_on * (t(increasing_phase) - t_start));
% Decreasing phase: Exponential decay from Kf1Max to Kf1Min
kf_end = Kf1Max - (Kf1Max - Kf1Min) * exp(Kf1Tau_on * (t_end - t_start));
decreasing_phase = t >= t_end;
kf_1(decreasing_phase) = Kf1Min + (kf_end - Kf1Min) * exp(Kf1Tau_off * (t(decreasing_phase) - t_end));
end
function kb_1 = calculate_kb(t, t_start, t_end, Kb1Max, Kb1Min, Kb1Tau_on, Kb1Tau_off)
% Initialize kb_1
kb_1 = zeros(size(t));
% Initial phase: Constant value Kb1Min
initial_phase = t < t_start;
kb_1(initial_phase) = Kb1Min;
% Increasing phase: Exponential growth from Kb1Min to Kb1Max
increasing_phase = t >= t_start & t < t_end;
kb_1(increasing_phase) = Kb1Max - (Kb1Max - Kb1Min) * exp(Kb1Tau_on * (t(increasing_phase) - t_start));
% Decreasing phase: Exponential decay from Kb1Max to Kb1Min
kb_end = Kb1Max - (Kb1Max - Kb1Min) * exp(Kb1Tau_on * (t_end - t_start));
decreasing_phase = t >= t_end;
kb_1(decreasing_phase) = Kb1Min + (kb_end - Kb1Min) * exp(Kb1Tau_off * (t(decreasing_phase) - t_end));
end
function dydt = ode_LR(t, y, kf_1, kb_1)
% Unpack the concentrations
Receptor_concentration = y(1);
C_LigandReceptor = y(2);
% Define the ODE system
dReceptor_dt = -kf_1(t) * Receptor_concentration + kb_1(t) * C_LigandReceptor;
d_CLigandReceptor_dt = kf_1(t) * Receptor_concentration - kb_1(t) * C_LigandReceptor;
% Pack the derivatives into a column vector
dydt = [dReceptor_dt; d_CLigandReceptor_dt];
end
0 件のコメント
回答 (1 件)
Harald
2024 年 1 月 22 日
Hi,
after running your code, combined_table does not contain anything called Kf1Max or even remotely similar to it.
Best wishes,
Harald
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1596906/image.png)
3 件のコメント
Image Analyst
2024 年 1 月 23 日
OK, fine - the values are already there in the table. It looks like Kf1Max is not even involved. Do you have an (x,y) coordinate for each value of R_peak? If so, then use surf
Harald
2024 年 1 月 23 日
@Ehtisham, it is of course always good to know where you ultimately want to go. However, I don't think this resolves the imminent problem: you are trying to extract a variable from a table that isn't there.
So, as Image Analyst suggests, you will need to obtain that variable from "somewhere" else.
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!