Facing difficulty to resolve code problem

% Define symbolic variables

syms ABA ABA_Receptor RAF_SnRK2_Cascade SnRK2_6 SnRK2_6_P B2_RAF B3_RAF t

syms k1 k2 k3 k4;

% Define kinetic parameters

k1 = 1; % Rate constant for ABA binding to receptor

k2 = 0.5; % Rate constant for RAF-SnRK2 interaction

k3 = 0.8; % Rate constant for SnRK2 dephosphorylation

k4 = 0.6; % Rate constant for SnRK2 transphosphorylation by RAFs

B2_RAF = 0.1; % Replace with actual values

B3_RAF = 0.2; % Replace with actual values

% Define the system of ODEs

eq1 = k1 * ABA - k2 * ABA_Receptor * RAF_SnRK2_Cascade;

eq2 = k2 * ABA_Receptor * RAF_SnRK2_Cascade - k3 * RAF_SnRK2_Cascade;

eq3 = k3 * RAF_SnRK2_Cascade - k4 * SnRK2_6 * (B2_RAF + B3_RAF);

eq4 = k4 * SnRK2_6 * (B2_RAF + B3_RAF);

% Convert the symbolic equations to a numerical function

odeSystem = matlabFunction([eq1; eq2; eq3; eq4], 'Vars', [t, ABA, ABA_Receptor, RAF_SnRK2_Cascade, SnRK2_6, SnRK2_6_P, B2_RAF, B3_RAF]);

% Initial conditions for the ODE solver

initial_conditions = [0, 0, 0, 0];

% Time span for simulation

tspan = 0:0.1:10;

% Solve the ODE numerically

[t, y] = ode45(@(t, y) odeSystem(t, y, ABA, ABA_Receptor, RAF_SnRK2_Cascade, SnRK2_6, SnRK2_6_P, B2_RAF, B3_RAF), tspan, initial_conditions);

% Extract numerical solutions

ABA_Receptor_sol = Y(:, 1);

RAF_SnRK2_Cascade_sol = Y(:, 2);

SnRK2_6_sol = Y(:, 3);

SnRK2_6_P_sol = Y(:, 4);

% Plot the results

figure;

subplot(2, 2, 1);

plot(t, ABA_Receptor_sol);

xlabel('Time');

ylabel('ABA Receptor');

subplot(2, 2, 2);

plot(t, RAF_SnRK2_Cascade_sol);

xlabel('Time');

ylabel('RAF-SnRK2 Cascade');

subplot(2, 2, 3);

plot(t, SnRK2_6_sol);

xlabel('Time');

ylabel('SnRK2.6');

subplot(2, 2, 4);

plot(t, SnRK2_6_P_sol);

xlabel('Time');

ylabel('Phosphorylated SnRK2.6');

subplot(2, 2, 2);

plot(t, RAF_SnRK2_Cascade_sol);

xlabel('Time');

ylabel('RAF-SnRK2 Cascade');

subplot(2, 2, 3);

plot(t, SnRK2_6_sol);

xlabel('Time');

ylabel('SnRK2.6');

subplot(2, 2, 4);

plot(t, SnRK2_6_P_sol);

xlabel('Time');

ylabel('Phosphorylated SnRK2.6');

Error using sym/matlabFunction>checkVars

Variable names must be valid MATLAB variable names.

Error in sym/matlabFunction (line 158)

vars = checkVars(funvars,opts);

Error in aba_SnRK (line 20)

odeSystem = matlabFunction([eq1; eq2; eq3; eq4], 'Vars', [t, ABA, ABA_Receptor, RAF_SnRK2_Cascade, SnRK2_6, SnRK2_6_P, B2_RAF, B3_RAF]);

[Ramtej]

Hi Shakeel,

The error you are encountering is due to the values being assigned to the input arguments "B2_RAF","B3_RAF" for your "odeSystem".

Do not assign the value to "B2_RAF", "B3_RAF" before creating the "odeFunction" using "matlabFunction".

[Dyuman Joshi]

Keeping the errors at one side, the 4 equations you have written do not depend on t or y.

How can you solve for variables that don't exist in the equations?

And in the same line, the ODE call does not make sense.

[Shakeel Ahmed]

Thanks for help

Answer 1

[Sam Chak]

Hi @Shakeel Ahmed

I recommend using the numerical approach directly when the state equations are available, as in your case. It is unnecessary to describe the ODEs in symbolic form and then convert to numerical form for the ode45() solver. For your problem, I mainly rearranged your code in a structure where it can be solved by the ode45() solver. All naming of the variables, parameters, and the system remains unchanged.

However, you didn't supply the value for a parameter named 'ABA', so please edit the value accordingly.

% Initial conditions for the ODE solver
initial_conditions = [0, 0, 0, 0];

% Time span for simulation
tspan = 0:0.1:10;

% Solve the ODE numerically
[t, y] = ode45(@odeSystem, tspan, initial_conditions);

% Extract numerical solutions
ABA_Receptor_sol = y(:, 1);
RAF_SnRK2_Cascade_sol = y(:, 2);
SnRK2_6_sol = y(:, 3);
SnRK2_6_P_sol = y(:, 4);

% Plot the results
figure;
subplot(2, 2, 1);
plot(t, ABA_Receptor_sol);
xlabel('Time');
ylabel('ABA Receptor');

subplot(2, 2, 2);
plot(t, RAF_SnRK2_Cascade_sol);
xlabel('Time');
ylabel('RAF-SnRK2 Cascade');

subplot(2, 2, 3);
plot(t, SnRK2_6_sol);
xlabel('Time');
ylabel('SnRK2.6');

subplot(2, 2, 4);
plot(t, SnRK2_6_P_sol);
xlabel('Time');
ylabel('Phosphorylated SnRK2.6');

function dydt = odeSystem(t, y)
    % Initialization
    dydt    = zeros(4, 1);

    % Definitions of the states
    ABA_Receptor        = y(1);
    RAF_SnRK2_Cascade   = y(2);
    SnRK2_6             = y(3);
    SnRK2_6_P           = y(4);

    % Define the parameters
    k1      = 1;    % Rate constant for ABA binding to receptor
    k2      = 0.5;  % Rate constant for RAF-SnRK2 interaction
    k3      = 0.8;  % Rate constant for SnRK2 dephosphorylation
    k4      = 0.6;  % Rate constant for SnRK2 transphosphorylation by RAFs
    B2_RAF  = 0.1;  % Replace with actual values
    B3_RAF  = 0.2;  % Replace with actual values
    ABA     = 1;    % this value is not provided by the OP
    
    % Define the system of ODEs
    dydt(1) = k1 * ABA - k2 * ABA_Receptor * RAF_SnRK2_Cascade;
    dydt(2) = k2 * ABA_Receptor * RAF_SnRK2_Cascade - k3 * RAF_SnRK2_Cascade;
    dydt(3) = k3 * RAF_SnRK2_Cascade - k4 * SnRK2_6 * (B2_RAF + B3_RAF);
    dydt(4) = k4 * SnRK2_6 * (B2_RAF + B3_RAF);
end