How can I fix ODE45 returning a straight line when using specific `tspan` points, instead of showing dynamic behavior over the time range?

How can I fix ODE45 returning a straight line when using specific `tspan` points, instead of showing dynamic behavior over the time range?

[Torsten]

By choosing more output points. This has already been excessively discussed with you.

[Torsten]

Use a loop:

Call ode45 from t=0 to t=500 with RLC = 0.1, restart ode45 from t=500 to t=1000 with the end values of the former integration and RLC = 0.5,...

[Torsten]

RLC = [ 0.1, 0.5, 10, 2];
Times = [ 0, 500, 1000, 2000];
y0 = ...;
T = [];
Y = [];
for i = 1:numel(Times)-1
  tspan = [Times(i) Times(i+1)];
  [t, y] = ode45(@(t,y)fun(t,y,RLC(i)),tspan,y0)
  T = [T;t];
  Y = [Y;y];
  y0 = y(end,:);
end

[Torsten]

The code below works, but I doubt it is really what you want. Try to explain yourself more clearly (maybe by using the Google Translator).

% Define parameters
Kf_Max = 100;  % maximum forward rate
Kb_Max = 80;  % maximum backward rate 
Kb_Min = 10 ;
RLC = [0.1, 0.5, 10, 1]; % RLC values
TauKf_ON = -0.1; % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF
TauKb_ON = -0.01;  % TauKb_ON
TauKb_OFF = -0.01; % TauKb_ON
PhaseTimes = [0,10, 1000, 1500, 2000];% phase times (each row defines a phase)

% 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)
[t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase (Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...
    Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, PhaseTimes,Non_Active_Receptor_concentration, Active_Receptor_concentration);

function [t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase(Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...
    Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, PhaseTimes,  Non_Active_Receptor_concentration, Active_Receptor_concentration)

% 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.01, 'RelTol', 1e-6, 'AbsTol', 1e-8);
t = [];
y = [];
for i = 1:numel(PhaseTimes)-1
    PhaseTime = [PhaseTimes(i) PhaseTimes(i+1)];
    % Define functions for forward and backward reaction rates
    Kf_L = @(t) calculate_kf(t, PhaseTime, Kf_Max, RLC(i), TauKf_ON, TauKf_OFF);
    Kb = @(t) calculate_kb(t, PhaseTime, Kb_Max, Kb_Min, RLC(i), TauKb_ON, TauKb_OFF);
    [T,Y] = ode45(@(t,y)ode_LR(t,y,Kf_L,Kb),PhaseTime,initial_conditions, options);
    initial_conditions = Y(end,:);
    t = [t;T];
    y = [y;Y];
end
% 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, 'b', '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 using element-wise division
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
% Nested function for calculating Kf_L based on time
function Kf_L = calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)
% Calculate maximum Kf_L based on RLC values using element-wise division
Kf_LMax = Kf_Max * (RLC ./ (1 + RLC));

% Default Kf_L to the maximum value of the first phase
Kf_L = Kf_LMax(1);

% Number of phases
num_phases = numel(RLC);

% Iterate through each phase to determine the current phase based on time
for j = 1:num_phases
    T_start = PhaseTimes(j); % Start time of the current phase
    if j < num_phases
        T_end = PhaseTimes(j + 1); % End time of the current phase
    else
        T_end = inf; % Handle last phase separately
    end
    
    % Check if the current time t is within the current phase
    if t >= T_start && t < T_end
        if j == 1
            % For the first phase, use the maximum value directly
            Kf_L = Kf_LMax(j);
        else
            % Time at the end of the previous phase
            prev_end = PhaseTimes(j-1);
            
            % Determine if RLC is increasing or decreasing and calculate Kf_L
            if RLC(j) > RLC(j - 1)
                % Increasing phase
                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_OFF * (t - prev_end));
            elseif  RLC(j) < RLC(j - 1)
                % Decreasing phase
                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_OFF * (t - prev_end));
            end
        end
        return; % Exit the function after finding the correct phase
    end
end
end
% Nested function for calculating Kb based on time
function Kb = calculate_kb(t, PhaseTimes, Kb_Max, Kb_Min, RLC, TauKb_ON, TauKb_OFF)
% Default to Kb_Max
Kb = Kb_Max;

% If all RLC values are the same, set Kb to Kb_Max and return
if all(RLC == RLC(1))
    Kb = Kb_Max;
    return;
end

% Iterate through each phase based on time
for j = 1:numel(RLC)
    T_start = PhaseTimes(j);   % Start time of the current phase
    if j < numel(PhaseTimes)
        T_end = PhaseTimes(j + 1); % End time of the current phase
    else
        T_end = inf; % For the last phase
    end
    
    % Check if current time is within the current phase
    if t >= T_start && t < T_end
        if j == 1
            % In the first phase, set Kb to Kb_Max
            Kb = Kb_Max;
        else
            prev_end = PhaseTimes(j-1); % Time at the end of the previous phase
            
            % Handle increasing or decreasing RLC values
            if RLC(j) > RLC(j - 1)
                % Increasing phase
                Kb_end = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
                Kb = Kb_Min + (Kb_end - Kb_Min) * exp(TauKb_OFF * (t - prev_end));
            elseif RLC(j) < RLC(j - 1)
                % Decreasing phase
                Kb_end = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));
                Kb = Kb_Min + (Kb_end - Kb_Min) * exp(TauKb_OFF * (t - prev_end));
            end
        end
        return; % Exit once the correct phase is found
    end
end
end
% ODE for calculation of the derivative 
function dydt = ode_LR(t, y, Kf_L, Kb)
Non_Active_Receptor_concentration = y(1);
Active_Receptor_concentration = y(2);

% Call the function handles for Kf_L and Kb
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
% Results of all the Phase Values 
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
% Customize function for Calculating the Peak values and time
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

[Torsten]

You ask me questions about the code that you made. How should I know better than you do ?

I set the initial condition for the new phase to the end value of the last phase. Maybe that's the problem.

for i = 1:numel(PhaseTimes)-1

PhaseTime = [PhaseTimes(i) PhaseTimes(i+1)];

% Define functions for forward and backward reaction rates

Kf_L = @(t) calculate_kf(t, PhaseTime, Kf_Max, RLC(i), TauKf_ON, TauKf_OFF);

Kb = @(t) calculate_kb(t, PhaseTime, Kb_Max, Kb_Min, RLC(i), TauKb_ON, TauKb_OFF);

[T,Y] = ode45(@(t,y)ode_LR(t,y,Kf_L,Kb),PhaseTime,initial_conditions, options);

initial_conditions = Y(end,:);

t = [t;T];

y = [y;Y];

end

Answer 1

[Steven Lord]

Your ode_LR function is defined to accept four input arguments, the two required arguments plus two additional parameters.

function dydt = ode_LR(t, y, Kf_L, Kb)

How are you using it, based on what's displayed in the error message?

Error in SimfileNPhase>@(t,y)ode_LR(t,y,RLC(i)) (line 20)
  [t, y] = ode45(@(t,y)ode_LR(t,y,RLC(i)),PhaseTime,initial_conditions, options);

Not only are you not calling it with enough input arguments (the anonymous function calls it with three inputs) the comment inside ode_LR implies that the Kf_L and Kb input arguments are function handles.

% Call the function handles for Kf_L and Kb
    dNonActiveReceptor_dt = -Kf_L(t) * Non_Active_Receptor_concentration + Kb(t) * Active_Receptor_concentration;

Is RLC(i), the value that you pass into ode_LR as the third input argument, a function handle? No, it is not. It's a numeric array. And t (the value that ode45 passes into your function as the first required input) is not guaranteed to be a valid index into a numeric array.