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

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]);

Answer 1

[Torsten]

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]);

[Ehtisham]

I tried but its giving me same error again.

[Torsten]

I ran the code with these two changes, and this error was healed.

There were NaN values in the results of the double loop - but that's a different story:

R*peak: 79.0231
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

[Torsten]

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]

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.