Clear Filters
Clear Filters

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

Answers (2)

Torsten
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.
  2 Comments
Ehtisham
Ehtisham on 4 Sep 2024
@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
Torsten on 4 Sep 2024
Edited: Torsten on 4 Sep 2024
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

Sign in to comment.


Rahul
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))
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>
  1 Comment
Ehtisham
Ehtisham on 5 Sep 2024
@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

Sign in to comment.

Tags

Products


Release

R2023a

Community Treasure Hunt

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

Start Hunting!