When i remove the element wise divison of the RLC beacuse RLC is single array it gives me this error Index exceeds the number of array elements. Index must not exceed 1.

When i remove the element wise divison of the RLC beacuse RLC is single array it gives me this error Index exceeds the number of array elements. Index must not exceed 1.

Need your guidance

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 10, 1];  % RLC values
TauKf_ON = -0.1; % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF
Kb_Max = 50;  % maximum backward rate  
PhaseTimes = [0, 100, 200, 300, 400];% phase times (each row defines a phase)
timespan = [0,  400];  % 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)

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

    % 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-4, 'AbsTol', 1e-6);

    % 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 / (1+ RLC));

    % 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, TauKf_OFF)
    % Calculate the maximum Kf_L for each phase based on RLC values
    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 based on the length of RLC
    num_phases = numel(RLC);

    % Iterate through each phase to determine the current phase based on time t
    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);

                % For phases other than the last one
                if j < num_phases
                    % Peak condition: RLC value is higher in the current phase compared to its neighbors
                    if 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_end - T_start));
                        Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j - 1)) * exp(TauKf_OFF * (t - prev_end));

                    % Increasing RLC condition: 
                    elseif RLC(j - 1) < RLC(j) && RLC(j) <= RLC(j + 1)
                        Kf_L = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_ON * (T_end - T_start));
                       

                    % Decreasing RLC condition: RLC value is lower in the current phase compared to neighbors
                    elseif RLC(j - 1) > RLC(j) && RLC(j) < RLC(j + 1)
                        Kf_L = Kf_LMax(j) + (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_OFF * (t - prev_end));

                    % Flat or decreasing RLC condition
                    elseif RLC(j - 1) > RLC(j) && RLC(j) >= RLC(j + 1)
                        Kf_L = Kf_LMax(j) + (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_OFF * (t - prev_end));
                    end
                else
                    % For the last phase, consider only the previous phase's value and time
                    if RLC(j - 1) < RLC(j)
                        Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_ON * (T_end - T_start));
                        Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j - 1)) * exp(TauKf_OFF * (t - prev_end));
                    elseif RLC(j - 1) > RLC(j)
                        Kf_L = Kf_LMax(j) + (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_ON * (t - prev_end));
                    end
                end
            end
            return;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kb = calculate_kb( Kb_Max)
    Kb = Kb_Max;  % Default to the maximum value
end


% function Kb = calculate_kb(t, PhaseTimes, Kb_Max, RLC)
%     % Default value for Kb is set to Kb_Max
%     Kb = Kb_Max; 
% 
%     % Check if all RLC values are the same
%     if all(RLC == RLC(1))
%         % If all RLC values are the same, set Kb to Kb_Max
%         Kb = Kb_Max; 
%         return;
%     end
% 
%     % Loop through each phase interval
%     for j = 1:numel(RLC)
%         % Define the start and end times for the current phase
%         T_start = PhaseTimes(j);
%         T_end = PhaseTimes(j + 1);
% 
%         % Check if current time falls within the current phase interval
%         if t >= T_start && t < T_end
%             % If it's the first phase, Kb is always Kb_Max
%             if j == 1
%                 Kb = Kb_Max;
%             else
%                 % Check if we're not at the last phase
%                 if j < numel(RLC)
%                     % Determine Kb based on the relative values of RLC
%                     % Compare current RLC value with previous and next values
%                     if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
%                         % Local maximum
%                         Kb = Kb_Max;
%                     elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
%                         % Local peak or plateau
%                         Kb = Kb_Max;
%                     elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
%                         % Local minimum
%                         Kb = Kb_Max;
%                     elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
%                         % Local dip or plateau
%                         Kb = Kb_Max;
%                     end
%                 else
%                     % Handle the case when j is the last phase
%                     if RLC(j-1) < RLC(j)
%                         % If last value is higher than previous
%                         Kb = Kb_Max;
%                     elseif RLC(j-1) > RLC(j)
%                         % If last value is lower than previous
%                         Kb = Kb_Max;
%                     end
%                 end
%             end
%             return; % Exit function once the appropriate Kb value is set
%         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;

% print the Active and NoN Active Receptor concentration 
fprintf('t = %f, Non_Active = %f, Active = %f, dNonActive/dt = %f, dActive/dt = %f\n', t, Non_Active_Receptor_concentration, Active_Receptor_concentration, dNonActiveReceptor_dt, dActiveReceptor_dt);




    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(Phase_time, Phase_concentration)
    % Compute the derivative 
    dt = diff(Phase_time);
    dy = diff(Phase_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(Phase_concentration(max_indices));
                % Find the lowest minimum
                [Rpeak_min, minIdx] = min(Phase_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 = Phase_time(max_indices(maxIdx));
                    peak_type = 'High';
                else
                    Rpeak = Rpeak_min;
                    Tpeak = Phase_time(min_indices(minIdx));
                    peak_type = 'Low';
                end
            % If only maxima exist
            elseif ~isempty(max_indices)
                [Rpeak, maxIdx] = max(Phase_concentration(max_indices));
                Tpeak = Phase_time(max_indices(maxIdx));
                peak_type = 'High';
            % If only minima exist
            elseif ~isempty(min_indices)
                [Rpeak, minIdx] = min(Phase_concentration(min_indices));
                Tpeak = Phase_time(min_indices(minIdx));
                peak_type = 'Low';
            end
        end
    end
end

Answer 1

[Torsten]

You must keep the elementwise division. So

Change

Kf_LMax = Kf_Max * (RLC / (1 + RLC));

to

Kf_LMax = Kf_Max * (RLC ./ (1 + RLC));

What else do you want Kf_LMax(i) to become than Kf_Max * RLC(i)/(1+RLC(i)) ?

And you should use ode15s instead of ode45 - it's way faster.

[Ehtisham]

@Torsten i want this Kf_LMax = Kf_Max * (RLC / (1 + RLC)); in code not

Kf_LMax = Kf_Max * (RLC ./ (1 + RLC)); i do not want the Element wise division so when i remove i get the error as i explained in question statment.

[Torsten]

Do you know what you get from

Kf_LMax = Kf_Max * (RLC / (1 + RLC));

if RLC is a vector ? You get nonsense (at least for your application).

RLC = [1 2 3 4 5];
Kf_Max = 1;
Kf_Max * (RLC / (1 + RLC))

ans = 0.7778

Kf_Max * sum(RLC.*(1+RLC))/sum((1+RLC).^2)

ans = 0.7778

Answer 2

[Rahul]

Hi Ehtisham,

I understand that you’re trying to solve a series of equations for a given RLC circuit and a set of initial conditions, where an error stating “Index exceeds the number of array elements. Index must not exceed 1” is thrown on program execution.

Specifically, in the function definition of ‘calculate_kft’, since variable ‘Kf_Max’ and division result involving vector ‘RLC’ are scalars themselves, the variable ‘Kf_LMax’ is being assigned a scalar value, which leads to the out of range index error statement:

RLC = [0.1, 0.5, 10, 1];
Kf_Max = 100;
x = (RLC / (1 + RLC))

x = 0.8786

Kf_LMax = Kf_Max * (RLC / (1 + RLC))

Kf_LMax = 87.8561

Further, variable Kb has been assigned a scalar value of ‘Kb_Max’ = 50, inside the function ‘calculate_kb’ definition. Later, these two variables have been assumed to be non-scalars, in the following functions:

• calculate_kf: variable ‘num_phases’ stores the size of ‘RLC’ vector, as 3. Hence, while inside the for-loop ‘Iterate through each phase to determine the current phase based on time ‘t’ using an index ‘j’, calculating ‘Kf_end’ and ‘Kf_L’ for ‘phases other than the last one’, ‘Kf_LMax’ is being indexed as a vector, causing the stated error.
• ode_LR: variable ‘Kb’ has been indexed using variable ‘t’, whereas earlier ‘Kb’ has been assigned a scalar value:

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;

You can either include element wise division while calculating Kf_LMax as shown below, or make appropriate modifications to the above-mentioned statements.

RLC = [0.1, 0.5, 10, 1];
Kf_Max = 100;
x = (RLC ./ (1 + RLC))

x = 1x4

0.0909 0.3333 0.9091 0.5000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Ehtisham]

@Rahul what kind of adjustment i need to ensure so i have a standard division in the forward reaction and also in the findpeak function. Kindly guide