rror using odearguments @(T,Y)ODE_​LR(T,Y,KF_​L,KB_1) returns a vector of length 2, but the length of initial conditions vector is 10002. The vector returned by @(T,Y)ODE_LR(T,Y

2 次查看（过去 30 天）
Ehtisham
Ehtisham 2024-4-18
编辑： Torsten 2024-4-19
Error using odearguments
@(T,Y)ODE_LR(T,Y,KF_L,KB_1) returns a vector of length 2, but the length of initial conditions vector is 10002. The vector returned
by @(T,Y)ODE_LR(T,Y,KF_L,KB_1) and the initial conditions vector must have the same number of elements.
Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in Corefile (line 48)
[t, y] = ode45(@(t, y) ode_LR(t, y, kf_L, kb_1), timespan, [Non_Active_Receptor_concentration; Active_Receptor_concentration]);

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2024-4-18
编辑：Torsten 2024-4-18
Use
Non_Active_Receptor_concentration0 = 100;
Active_Receptor_concentration0 = 0;
instead of
Non_Active_Receptor_concentration = 100;
Active_Receptor_concentration = 0;
and
[t, y] = ode45(@(t, y) ode_LR(t, y, kf_L, kb_1), timespan, [Non_Active_Receptor_concentration0; Active_Receptor_concentration0]);
instead of
[t, y] = ode45(@(t, y) ode_LR(t, y, kf_L, kb_1), timespan, [Non_Active_Receptor_concentration; Active_Receptor_concentration]);
  4 个评论
显示 2更早的评论
Torsten
Torsten 2024-4-18
编辑：Torsten 2024-4-18
Ran in:
I used ode15s because it's faster for MATLAB online, but it will work equally well with ode45.
As said - NaN values appear in the evaluation of the results, but that's a different story. I'm confident that you will be able to solve this problem on your own.
Corefile()
R*peak: 79.023 T-peak: 2.7 50% form peak to ss : NaN R* Steady state: NaN T-50 : -2.7 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 88.2817 T-peak: 2.7 50% form peak to ss : NaN R* Steady state: NaN T-50 : -2.7 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 91.8713 T-peak: 2.7 50% form peak to ss : NaN R* Steady state: NaN T-50 : -2.7 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 93.7763 T-peak: 2.7 50% form peak to ss : NaN R* Steady state: NaN T-50 : -2.7 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 94.9596 T-peak: 2.7 50% form peak to ss : NaN R* Steady state: NaN T-50 : -2.7 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 75.9269 T-peak: 2.1 50% form peak to ss : NaN R* Steady state: NaN T-50 : -2.1 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 86.3186 T-peak: 2.1 50% form peak to ss : NaN R* Steady state: NaN T-50 : -2.1 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 90.4439 T-peak: 2.1 50% form peak to ss : NaN R* Steady state: NaN T-50 : -2.1 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 92.6545 T-peak: 2.2 50% form peak to ss : NaN R* Steady state: NaN T-50 : -2.2 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 94.0367 T-peak: 2.1 50% form peak to ss : NaN R* Steady state: NaN T-50 : -2.1 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 73.4103 T-peak: 1.8 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.8 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 84.6695 T-peak: 1.8 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.8 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 89.2276 T-peak: 1.8 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.8 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 91.6979 T-peak: 1.8 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.8 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 93.2452 T-peak: 1.8 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.8 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 71.2471 T-peak: 1.7 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.7 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 83.2035 T-peak: 1.7 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.7 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 88.1387 T-peak: 1.6 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.6 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 90.833 T-peak: 1.6 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.6 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 92.529 T-peak: 1.6 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.6 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 69.3119 T-peak: 1.5 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.5 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 81.8749 T-peak: 1.5 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.5 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 87.1415 T-peak: 1.5 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.5 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 90.0354 T-peak: 1.5 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.5 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN R*peak: 91.8659 T-peak: 1.5 50% form peak to ss : NaN R* Steady state: NaN T-50 : -1.5 R*peak/R*ss: NaN Δ (R*peak-R*ss) : NaN
function Corefile()
timespan = 0:0.1:500;
Non_Active_Receptor_concentration0 = 100;
Active_Receptor_concentration0 = 0;
Kb1Max_values = [100, 200, 300, 400, 500];
Kf1Max_values = [100 , 200, 300, 400, 500];
% Create a cell array to store tables for different Kf1Max and Kb1Max ranges
tables_cell = cell(length(Kb1Max_values), length(Kf1Max_values));
% Initialize arrays to store Kf1Max and Kb1Max values for each iteration
Kb1Max_result = zeros(length(Kb1Max_values), length(Kf1Max_values));
Kf1Max_result = zeros(length(Kb1Max_values), length(Kf1Max_values));
R_peak_values = zeros(length(Kb1Max_values), length(Kf1Max_values));
T_peak_values = zeros(length(Kb1Max_values), length(Kf1Max_values));
R_steady_state_values = zeros(length(Kb1Max_values), length(Kf1Max_values));
T_50_values = zeros(length(Kb1Max_values), length(Kf1Max_values));
R_peak_R_ss_values = zeros(length(Kb1Max_values), length(Kf1Max_values));
Delta_values = zeros(length(Kb1Max_values), length(Kf1Max_values));
for i = 1:length(Kb1Max_values)
Kb1Max = Kb1Max_values(i);
for j = 1:length(Kf1Max_values)
Kf1Max = Kf1Max_values(j);
T_startLigand = 100;
T_endLigand = 1000;
TauKFON = -1;
TauKFOFF= -1;
Kb1Min = 10;
TauKBON = -0.01;
TauKBOFF = -0.01;
L_min = 0.01;
L_Active = 1 ;
% Define functions for forward and backward reaction rates while
% forward reaction depends on the Time and Ligand too
kf_L = @(t) calculate_kf(t, T_startLigand, T_endLigand, Kf1Max, L_Active, L_min, TauKFON, TauKFOFF);
kb_1 = @(t) calculate_kb(t, T_startLigand, T_endLigand, Kb1Max, Kb1Min, TauKBON, TauKBOFF);
% Initialize KF_LMaxA with the minimum amount of the Ligand
KF_LMaxA = Kf1Max* (L_min / (L_min + 1)); % KF_LMaxA is calculated from KfMax means (steady state with [L]= Lmin)
% Initialize KF_LMaxB
KF_LMaxB = Kf1Max* (L_Active / (L_Active + 1)); % KF_LMaxB is calculated from KfMax means (Activation with [L] = LActive)
[t, y] = ode15s(@(t, y) ode_LR(t, y, kf_L, kb_1), timespan, [Non_Active_Receptor_concentration0; Active_Receptor_concentration0]);
% Extract the concentrations
Non_Active_Receptor_concentration = y(:, 1);
Active_Receptor_concentration = y(:, 2);
% Calculating the R*peak (peak value means maximum concentration of active receptors, and peak_time_idx contains the index of the time at which this peak concentration occurs and max (Active receptor concentartion ) means getting maximum value in the vector).
[peak_value, peak_time_idx] = max(Active_Receptor_concentration(t >= T_startLigand & t <= T_endLigand));
% This line find function is used to locate the first index in the t vector that is within the specified time range. The -1 is added to convert the index from the subarray back to the original array, and the peak_time_idx is then added to get the correct index within the entire t vector.
peak_time_idx = find(t >= T_startLigand & t <= T_endLigand, 1, 'first') - 1 + peak_time_idx;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculating the T-Peak
% This line finds the time duration from the start of ligand stimulation to the time at which the active receptor concentration reaches its peak. This is done by subtracting the time at the start of ligand stimulation from the time at the peak.
time_to_peak = t(peak_time_idx)-T_startLigand;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the R*ss
% This line find the differnce last 10 points till the Tendligand
steady_state = mean(Active_Receptor_concentration( t>=(T_endLigand-10) & t<=T_endLigand));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calclulating the T-50 explain every step in the comment
% This line calculates the midpoint between the steady_state and the peak_value.
% It is essentially finding the value that is halfway between the baseline (steady state) and the peak value of the Active_Receptor_concentration.
half_peak_value = (steady_state + peak_value) / 2;
% This line finds the index (idx_50_percent) where the Active_Receptor_concentration is closest to the calculated half_peak_value.
% The min(abs(...)) part is used to find the minimum absolute difference between the Active_Receptor_concentration values and the half_peak_value.
[~, idx_50_percent] = min(abs(Active_Receptor_concentration - half_peak_value));
% This line calculates the time it takes for the receptor concentration to reach 50% of the peak value.
% It subtracts the time corresponding to the peak_time_idx from the time corresponding to the index idx_50_percent.
time_to_50_percent = t(idx_50_percent)-time_to_peak ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculating the ratio of the R*peak/R* steady state
peak_to_steady_state_ratio = peak_value / steady_state;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculating 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)]);
% Store Kf1Max and Kb1Max values for each iteration
Kf1Max_result(i, j) = Kf1Max;
Kb1Max_result(i, j) = Kb1Max;
R_peak_values(i, j) = peak_value;
T_peak_values(i, j) = time_to_peak;
R_steady_state_values(i, j) = steady_state;
T_50_values(i, j) = time_to_50_percent;
R_peak_R_ss_values(i, j) = peak_to_steady_state_ratio;
Delta_values(i, j) = Delta;
% Combine all results into a single table
combined_table = table(Kf1Max_result(:), Kb1Max_result(:), R_peak_values(:), T_peak_values(:), R_steady_state_values(:), T_50_values(:), R_peak_R_ss_values(:), Delta_values(:), ...
'VariableNames', {'Kb1Max', 'Kf1Max', 'R_peak', 'T_peak', 'R_steady_state', 'T_50', 'R_peak_R_ss', 'Delta'});
% Generate a unique filename for the combined table
combined_csv_filename = 'combined_interaction_data_all.csv';
writetable(combined_table, combined_csv_filename);
end
end
function kf_L = calculate_kf(t, T_startLigand, T_endLigand, Kf1Max, L_Active, L_min, TauKFON, TauKFOFF)
% Initial phase: Constant value KF_LMaxA
if (t < T_startLigand)
kf_L = KF_LMaxA;
% Increasing phase: Exponential growth from KF_LMaxA to KF_LMaxB
elseif (t >= T_startLigand) && (t < T_endLigand)
kf_L = KF_LMaxB - (KF_LMaxB - KF_LMaxA) * exp(TauKFON * (t - T_startLigand));
% Decreasing phase: Exponential decay from KF_LMaxA to KF_LMaxB
else
kf_end = KF_LMaxB - (KF_LMaxB - KF_LMaxA) * exp(TauKFON * (T_endLigand - T_startLigand));
kf_L = KF_LMaxA + (kf_end - KF_LMaxA) * exp(TauKFOFF * (t - T_endLigand));
end
end
function kb_1 = calculate_kb(t, T_startLigand, T_endLigand, Kb1Max, Kb1Min, TauKBON, TauKBOFF)
% Initial phase: Constant value Kb1Min
if (t < T_startLigand)
kb_1 = Kb1Min;
% Increasing phase: Exponential growth from Kb1Min to Kb1Max
elseif (t >= T_startLigand) && (t < T_endLigand)
kb_1 = Kb1Max - (Kb1Max - Kb1Min) * exp(TauKBON * (t - T_startLigand));
% Decreasing phase: Exponential decay from Kb1Max to Kb1Min
else
kb_end = Kb1Max - (Kb1Max - Kb1Min) * exp(TauKBON * (T_endLigand - T_startLigand));
kb_1 = Kb1Min + (kb_end - Kb1Min) * exp(TauKBOFF * (t - T_endLigand));
end
end
function dydt = ode_LR(t, y, kf_L, kb_1)
% Unpack the concentrations
Non_Active_Receptor_concentration = y(1);
Active_Receptor_concentration = y(2);
% Define the ODE system
dNonActiveReceptor_dt = -kf_L(t) * Non_Active_Receptor_concentration + kb_1(t) * Active_Receptor_concentration;
dActiveReceptor_dt = kf_L(t) * Non_Active_Receptor_concentration - kb_1(t) * Active_Receptor_concentration;
% Pack the derivatives into a column vector
dydt = [dNonActiveReceptor_dt; dActiveReceptor_dt];
end
% Generate 3D mesh plot for R* Peak Values
figure;
surf(Kf1Max_result, Kb1Max_result, R_peak_values);
xlabel('Kf1Max');
ylabel('Kb1Max');
zlabel('R* Peak Value');
title(' R* Peak Values');
grid on;
colorbar;
% Set colormap to 'jet' for better visualization
colormap('jet');
% Set z-axis range for the plot
zlim([0, 100]);
% Generate 3D mesh plot for R* Steady State Values
figure;
surf(Kf1Max_result, Kb1Max_result, R_steady_state_values);
xlabel('Kf1Max');
ylabel('Kb1Max');
zlabel('R* Steady State Value');
title(' R* Steady State Values');
grid on;
colorbar;
% Set colormap to 'jet' for better visualization
colormap('jet');
% Set z-axis range for the plot
zlim([0, 100]);
% Generate 3D mesh plot for T-50 Values
figure;
surf(Kf1Max_result, Kb1Max_result, T_50_values);
xlabel('Kf1Max');
ylabel('Kb1Max');
zlabel('T-50 Value');
title(' T-50 Values');
grid on;
colorbar;
% Set colormap to 'jet' for better visualization
colormap('jet');
% Set z-axis range for the plot
zlim([0, 100]);
% Generate 3D mesh plot for T-peak Values
figure;
surf(Kf1Max_result, Kb1Max_result, T_peak_values);
xlabel('Kf1Max');
ylabel('Kb1Max');
zlabel('T-peak Value');
title(' T-peak Values');
grid on;
colorbar;
% Set colormap to 'jet' for better visualization
colormap('jet');
% Set z-axis range for the plot
zlim([0, 5]);
% Generate 3D mesh plot for R*peak/R*ss Values
figure;
surf(Kf1Max_result, Kb1Max_result, R_peak_R_ss_values);
xlabel('Kf1Max');
ylabel('Kb1Max');
zlabel('R*peak/R*ss Value');
title(' R*peak/R*ss Values');
grid on;
colorbar;
% Set colormap to 'jet' for better visualization
colormap('jet');
% Set z-axis range for the plot
zlim([0, 5]);
% Generate 3D mesh plot for Delta Values
figure;
surf(Kf1Max_result, Kb1Max_result, Delta_values);
xlabel('Kf1Max');
ylabel('Kb1Max');
zlabel('Δ = R*peak - R*ss');
title(' Δ = R*peak - R*ss ');
grid on;
colorbar;
% Set z-axis range for the plot
zlim([0, 100]);
% Set colormap to 'jet' for better visualization
colormap('jet');
end
Torsten
Torsten 2024-4-19
编辑：Torsten 2024-4-19
Ran in:
Active_Receptor_concentration( t>=(T_endLigand-10) & t<=T_endLigand)
is empty, and
mean([])
ans = NaN
gives NaN. This explains steady_state = NaN.
Most probably, the other quantities follow.
Your timespan must be at least up to T_endLigand.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by