# 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.

10 views (last 30 days)

Show older comments

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

##### 0 Comments

### Answers (2)

Torsten
on 4 Sep 2024

Edited: Torsten
on 4 Sep 2024

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.

Rahul
on 5 Sep 2024

Edited: Rahul
on 5 Sep 2024

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

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

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

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!