Selecting appropriate values for the parameters in the Gierer-Meinhardt activator-inhibitor model

For the following activator-inhibitor system (Alfred Gierer-Hans Meinhardt model),

with all constants positive and initial conditions , examine the effect of the parameters () on the behavior of concentrations . Initially, we have to consider .

My problem is that I don't undestand how to choose the right parameters every time. For example but the code works and for or but when I chose you can see below, the results are not correct.

Could anyone explain please? Thanks in advance

• Parameter k : Range: 0.01 to 1

% Define common parameters
    P = 0.75;   % Production rate of u
    a = 1;  % rate of u inhibition by v
    b = 0.9; % Production rate of v stimulated by u
    mu = 0.5; % Decay rate of u
    ks = [0.01, 0.1, 1];
    T = 10;  % Total simulation time
    dt = 0.01; % Time step
    t = 0:dt:T; % Time vector

        % Initial conditions
    u0 = 1; % Initial concentration of u
    v0 = 0.5; % Initial concentration of v
%% Initialize a figure
figure(1);
hold on

%% Loop through each k value
for k = ks
    N = length(t); % Number of time steps
    % Initialize arrays for u and v
    u = zeros(1, N);
    v = zeros(1, N);
    % Set initial conditions
    u(1) = u0;
    v(1) = v0;
    % Simulate the system's response
    for i = 2:N
        % Calculate the changes in u and v using the system equations
        f1 = P - a * (u(i-1)^2 / v(i-1)) - mu * u(i-1); % Change in activator
        f2 = b * u(i-1)^2 - k * v(i-1); % Change in inhibitor
        % Update u and v using Euler's method
        u(i) = u(i-1) + dt * f1;
        v(i) = v(i-1) + dt * f2;
    end
    % Plot the concentration of u and v for the current k value
    plot(t, u, 'DisplayName', ['u(t), k = ', num2str(k)]);
    plot(t, v, '--', 'DisplayName', ['v(t), k = ', num2str(k)]);
end

%% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Concentration');
title('Dynamics of u and v over Time for Different k Values');
legend('show');
grid on;
hold off;

• Parameter μ : Range: 0.01 to 1

% Define common parameters
P = 0.75;   % Production rate of u
a = 1;  % Rate of u inhibition by v
b = 0.9; % Production rate of v stimulated by u
k = 1;  % Decay rate of v (fixed in this scenario)
mus = [0.01, 0.1, 9];  % Different decay rates of u to test
T = 10;  % Total simulation time
dt = 0.01; % Time step
t = 0:dt:T; % Time vector

% Initial conditions
u0 = 1; % Initial concentration of u
v0 = 0.5; % Initial concentration of v

%% Initialize a figure
figure(2);
hold on

%% Loop through each mu value
for mu = mus
    N = length(t); % Number of time steps
    % Initialize arrays for u and v
    u = zeros(1, N);
    v = zeros(1, N);
    % Set initial conditions
    u(1) = u0;
    v(1) = v0;
    % Simulate the system's response
    for i = 2:N
        % Calculate the changes in u and v using the system equations
        f1 = P - a * (u(i-1)^2 / v(i-1)) - mu * u(i-1); % Change in activator
        f2 = b * u(i-1)^2 - k * v(i-1); % Change in inhibitor
        % Update u and v using Euler's method
        u(i) = u(i-1) + dt * f1;
        v(i) = v(i-1) + dt * f2;
    end
    % Plot the concentration of u and v for the current mu value
    plot(t, u, 'DisplayName', ['u(t), mu = ', num2str(mu)]);
    plot(t, v, '--', 'DisplayName', ['v(t), mu = ', num2str(mu)]);
end

%% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Concentration');
title('Dynamics of u and v over Time for Different mu Values');
legend('show');
grid on;
hold off;

• Parameter b : Range: 0.01 to 10

   
   % Define common parameters
P = 0.75;   % Production rate of u
a = 1;  % Rate of u inhibition by v
mu = 0.5; % Decay rate of u
k = 1;  % Decay rate of v 
bs = [0.1, 1, 9.9];  % Different rates for b to test
T = 10;  % Total simulation time
dt = 0.01; % Time step
t = 0:dt:T; % Time vector

% Initial conditions
u0 = 1; % Initial concentration of u
v0 = 0.5; % Initial concentration of v

%% Initialize a figure
figure(3);
hold on

%% Loop through each b value
for b = bs
    N = length(t); % Number of time steps
    % Initialize arrays for u and v
    u = zeros(1, N);
    v = zeros(1, N);
    % Set initial conditions
    u(1) = u0;
    v(1) = v0;
    % Simulate the system's response
    for i = 2:N
        % Calculate the changes in u and v using the system equations
        f1 = P - a * (u(i-1)^2 / v(i-1)) - mu * u(i-1); % Change in activator
        f2 = b * u(i-1)^2 - k * v(i-1); % Change in inhibitor
        % Update u and v using Euler's method
        u(i) = u(i-1) + dt * f1;
        v(i) = v(i-1) + dt * f2;
    end
    % Plot the concentration of u and v for the current b value
    plot(t, u, 'DisplayName', ['u(t), b = ', num2str(b)]);
    plot(t, v, '--', 'DisplayName', ['v(t), b = ', num2str(b)]);
end

%% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Concentration');
title('Dynamics of u and v over Time for Different b Values');
legend('show');
grid on;
hold off;

• Parameter a : Range: 0.1 to 10

% Define common parameters

P = 0.75;   % Production rate of u
b = 0.9; % Production rate of v stimulated by u
mu = 0.5; % Decay rate of u
k = 1;  % Decay rate of v (constant in this scenario)
as = [0.1, 1, 9.9];  % Different rates of a to test
T = 10;  % Total simulation time
dt = 0.01; % Time step
t = 0:dt:T; % Time vector

% Initial conditions
u0 = 1; % Initial concentration of u
v0 = 0.5; % Initial concentration of v

%% Initialize a figure
figure(4);
hold on

%% Loop through each a value
for a = as
    N = length(t); % Number of time steps
    % Initialize arrays for u and v
    u = zeros(1, N);
    v = zeros(1, N);
    % Set initial conditions
    u(1) = u0;
    v(1) = v0;
    % Simulate the system's response
    for i = 2:N
        % Calculate the changes in u and v using the system equations
        f1 = P - a * (u(i-1)^2 / v(i-1)) - mu * u(i-1); % Change in activator
        f2 = b * u(i-1)^2 - k * v(i-1); % Change in inhibitor
        % Update u and v using Euler's method
        u(i) = u(i-1) + dt * f1;
        v(i) = v(i-1) + dt * f2;
    end
    % Plot the concentration of u and v for the current a value
    plot(t, u, 'DisplayName', ['u(t), a = ', num2str(a)]);
    plot(t, v, '--', 'DisplayName', ['v(t), a = ', num2str(a)]);
end

%% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Concentration');
title('Dynamics of u and v over Time for Different a Values');
legend('show');
grid on;
hold off;

Answer 1

[Sam Chak]

Hi @Athanasios Paraskevopoulos

As usual, we begin by identifying the critical point(s) of the system, if they exist. However, I believe your code is correct, although it may be less accurate due to the use of the Euler method. If you examine my results generated by the ode45 solver (with ), you will notice that the state rapidly increases to a very high value before eventually converging to the calculated critical point. It's also worth noting that I selected different initial values to avoid encountering a division-by-zero singularity.

syms u v

%% parameters
P       = 0.75; % Production rate of u
b       = 0.9;  % Production rate of v stimulated by u
mu      = 0.5;  % Decay rate of u
k       = 1;    % Decay rate of v (constant in this scenario)
a       = 9.9;  % Different rates of a to test

%% find critical points
eq1     = P - a*u^2/v - mu*u == 0;
eq2     = b*u^2 - k*v == 0;
eqn     = [eq1, eq2];
sol     = solve(eqn);
ucrit   = double(sol.u)

ucrit = -20.5000

vcrit   = double(sol.v)

vcrit = 378.2250

%% call ode45 solver and plot result
tspan   = [0 100];
x0      = [1; 0.5];
options = odeset('RelTol', 1e-4, 'AbsTol', 1e-6);
[t, x]  = ode45(@(t, x) ode(t, x, P, b, mu, k, a), tspan, x0, options);
plot(t, x), grid on, legend('u', 'v'), xlabel('t')

%% Check if the states (u & v) settle near the critical point (equilibrium)
x(end,:)

ans = 1x2

-20.5002 378.2246

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

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

%% Gierer-Meinhardt activator-inhibitor model
function dx = ode(t, x, P, b, mu, k, a)
    u       = x(1);
    v       = x(2);
    dx(1,1) = P - a*u^2/v - mu*u;
    dx(2,1) = b*u^2 - k*v;
end

[Athanasios Paraskevopoulos]

Hello! Thank you for your answer!

When you said that you chose different initial values you mean this block of code

    options = odeset('RelTol', 1e-4, 'AbsTol', 1e-6);

[Sam Chak]

Nope. I meant the initial values are chosen at and , when .

x0      = [1; 0.5];

[Athanasios Paraskevopoulos]

Sorry but I didn’t get it. How are they different than mine?

[Sam Chak]

Previously, I didn't verify the initial values you selected. It appears that we both chose the same values for a = 9.9. However, I randomly selected those values to avoid encountering a singularity. I have a tendency to select x1(0) = 1.

[Athanasios Paraskevopoulos]

Hello again! Thank you for your explanation! I tried this to my code.

I change the time step to

Also, I used of `epsilon` to prevent division by zero when calculating f1

But as you can see our plot are different on the interval [0,100]. Does this happens because Euler's method?

% Define common parameters
P = 0.75;   % Production rate of u
b = 0.9; % Production rate of v stimulated by u
mu = 0.5; % Decay rate of u
k = 1;  % Decay rate of v (constant in this scenario)
as =  [9.9];  % Different rates of a to test
T = 100;  % Total simulation time
dt = 0.005; % Time step
t = 0:dt:T; % Time vector
epsilon = 1e-6; % Small number to prevent division by zero


% Initial conditions
u0 = 1; % Initial concentration of u
v0 = 0.5; % Initial concentration of v


%% Initialize a figure
figure(4);
hold on
%% Loop through each a value
for a = as
    N = length(t); % Number of time steps
    % Initialize arrays for u and v
    u = zeros(1, N);
    v = zeros(1, N);
    % Set initial conditions
    u(1) = u0;
    v(1) = v0;
    % Simulate the system's response
    for i = 2:N
        % Calculate the changes in u and v using the system equations
        f1 = P - a * (u(i-1)^2 / (v(i-1) + epsilon)) - mu * u(i-1); % Change in activator
        f2 = b * u(i-1)^2 - k * v(i-1); % Change in inhibitor
        % Update u and v using Euler's method
        u(i) = u(i-1) + dt * f1;
        v(i) = v(i-1) + dt * f2;
    end
    % Plot the concentration of u and v for the current a value
    plot(t, u, 'DisplayName', ['u(t), a = ', num2str(a)]);
    plot(t, v, '--', 'DisplayName', ['v(t), a = ', num2str(a)]);
end

%% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Concentration');
title('Dynamics of u and v over Time for Different a Values');
legend('show');
grid on;
hold off;

[Torsten]

If v becomes 0, the ODE has a singularity. You cannot just heal it by adding epsilon.

The results after the singularity will become almost arbitrary if you keep on integrating.

The same is true when applying ode45 as done above. The solver should stop when the singularity is reached. That this is not the case shows that also sophisticated solvers sometimes fail.

[Athanasios Paraskevopoulos]

If I understood correctly I can not do anything to fix this problem. Correct?

[Sam Chak]

Hi @Athanasios Paraskevopoulos

What is your objective in varying the parameters in the Gierer-Meinhardt activator-inhibitor model? By observing the model's behavior and drawing conclusions, one can make informed choices for initial values within a radius from the asymptotically stable region. Take a look at my example below with , .

It's crucial to recognize that the parameters impact the position of the critical point. Are you aiming to determine the five parameters {P, b, μ, k, a} such that the control-free trajectory settles at the desired critical point? What is your desired critical point?

Alternatively, if you have control over the model, you can identify a suitable location to introduce a manipulated input variable (ϵ) that allows the output variables (activator, u and inhibitor, v) to reach the desired operating point. This is commonly seen in PID-controlled drug delivery systems. However, injecting ϵ at the denominator, as proposed in your model, is not a wise choice.

syms u v

%% parameters
P       = 0.75; % Production rate of u
b       = 0.9;  % Production rate of v stimulated by u
mu      = 0.5;  % Decay rate of u
k       = 1;    % Decay rate of v (constant in this scenario)
a       = 9.9;  % Different rates of a to test

%% find critical points
eq1     = P - a*u^2/v - mu*u == 0;
eq2     = b*u^2 - k*v == 0;
eqn     = [eq1, eq2];
sol     = solve(eqn);
ucrit   = double(sol.u)

ucrit = -20.5000

vcrit   = double(sol.v)

vcrit = 378.2250

%% call ode45 solver and plot result
tspan   = [0 40];
x0      = [-10; 200];   % u0 = -10 and v0 = 200
options = odeset('RelTol', 1e-4, 'AbsTol', 1e-6);
[t, x]  = ode45(@(t, x) ode(t, x, P, b, mu, k, a), tspan, x0, options);
plot(t, x), grid on, legend('u', 'v'), xlabel('t')

%% Check if the states (u & v) settle near the critical point (equilibrium)
x(end,:)

ans = 1x2

-20.5009 378.2002

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

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

%% Gierer-Meinhardt activator-inhibitor model
function dx = ode(t, x, P, b, mu, k, a)
    u       = x(1);
    v       = x(2);
    dx(1,1) = P - a*u^2/v - mu*u;
    dx(2,1) = b*u^2 - k*v;
end

[Athanasios Paraskevopoulos]

Hello @Sam Chak! My aim is to see the effect of the parameters on the system for different values.

[Sam Chak]

Hi @Athanasios Paraskevopoulos,

Thank you for your clarification.

1. Have you discovered any definitive findings regarding the effect of the parameters using the Euler method?
2. If you switch to the Dormand–Prince method (ode45), do the findings change?
3. You mentioned about fixing the singularity problem with epsilon. How is singularity interpreted in your system? For example, what does it physically mean when the concentration of the inhibitor (v) is zero?

I'm unsure if the initial values are also important in your study. Considering the initial values and excluding P, there are a total of six parameters (a, b, μ, k, , ). If your study is limited to practical, non-negative concentrations of u and v, it would be advisable to calculate the critical point and select initial values within its neighborhood, so long as v doesn't decrease to 0.

[Torsten]

@Athanasios Paraskevopoulos

If I understood correctly I can not do anything to fix this problem. Correct?

The problem seems inherent to the ODE system. In such cases, the integrator usually stops near the singularity with the warning that the minimum stepsize is reached.

If you have an ODE with tan(x) as solution, e.g., you also cannot make the solver jump over pi/2.

[Sam Chak]

Hi @Torsten,

I'd like to inquire about cases involving removable singularities. In such situations, where the division-by-zero expression is undefined but has a limit as the denominator approaches 0, do we need to halt the ode45 solver and restart the integration when the state crosses zero?

dx      = @(t, x) - (sinh(x)/x - 1)*sign(x) + 0.05*sin(2*pi/20*t);
tspan   = [0 200];
x0      = 1;
[t, x]  = ode45(dx, tspan, x0);
plot(t, x), grid on, xlabel('t'), ylabel('x(t)')

[Torsten]

If x hits exactly 0 in your case from above, it does not help that you have a removable singularity or an essential one - you will get a division by zero.

The point is that the function (sinh(x)/x - 1) remains bounded near x = 0 - thus your function will be treated as if you were using

if x ~=0
    (sinh(x)/x - 1)
else
    0
end

because chances are very low that x = 0 is hit exactly during integration.

Different for tan(x):

tspan = [0 pi];
fun = @(t,y) 1/cos(t)^2;
[T,Y] = ode45(fun,tspan,0);

Warning: Failure at t=1.570796e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

plot(T,Y)

I guess in the case given chances are low that the singularity is removable since both u and v had to be 0 is the point of interest.

[Sam Chak]

Hi @Torsten

Thank you for your explanation. I will make sure to exercise caution when encountering various types of singularities.

[Athanasios Paraskevopoulos]

Hi@Sam Chak! My professor sau that we have to use Euler Method only for this project!

The exercise gives us that . And you have guessed correctly that the study is limited to practical, non-negative concentrations of u and v and we want to create a code to visualize the impact of parameter

[Athanasios Paraskevopoulos]

@Sam Chak @Torsten I think the following code is what I need! I have a block of code that insures that v does not decrease to zero or become negative

% Define common parameters
P = 0.75;   % Production rate of u
b = 0.9; % Production rate of v stimulated by u
mu = 0.5; % Decay rate of u
k = 1;  % Decay rate of v (constant in this scenario)
as = [0.1, 2, 4];  % Different rates of a to test
T = 10;  % Total simulation time
dt = 0.01; % Time step
t = 0:dt:T; % Time vector

% Initial conditions
u0 = 1; % Initial concentration of u
v0 = 0.5; % Initial concentration of v

%% Initialize a figure
figure(1);
hold on

%% Loop through each a value
for a = as
    N = length(t); % Number of time steps
    % Initialize arrays for u and v
    u = zeros(1, N);
    v = zeros(1, N);
    % Set initial conditions
    u(1) = u0;
    v(1) = v0;
    
    % Simulate the system's response
    for i = 2:N
        % Ensure v does not decrease to zero or become negative
        v_prev = max(v(i-1), 1e-5); % Prevent division by zero or negative v
        
        % Calculate the changes in u and v using the system equations
        f1 = P - a * (u(i-1)^2 / v_prev) - mu * u(i-1); % Change in activator
        f2 = b * u(i-1)^2 - k * v_prev; % Change in inhibitor
        
        % Update u and v using Euler's method
        u(i) = u(i-1) + dt * f1;
        v(i) = v(i-1) + dt * f2;
    end
    % Plot the concentration of u and v for the current a value
    plot(t, u, 'DisplayName', ['u(t), a = ', num2str(a)]);
    plot(t, v, '--', 'DisplayName', ['v(t), a = ', num2str(a)]);
end

%% Add labels, legend, and grid to the plot
xlabel('Time (s)');
ylabel('Concentration');
title('Dynamics of u and v over Time for Different a Values');
legend('show');
grid on;
hold off;

[Torsten]

As said, artificially keeping v > 0 will give you wrong results for the ODE if in reality, v would become 0. But in the time span up to 10 sec, setting

f1 = P - a * (u(i-1)^2 / v(i-1)) - mu * u(i-1); % Change in activator

instead of

f1 = P - a * (u(i-1)^2 / v_prev) - mu * u(i-1); % Change in activator

doesn't seem to change the results.

[Athanasios Paraskevopoulos]

@Torsten So what should I do with this method?

[Torsten]

What "method" do you mean ? If a differential equation has a singularity at some t-value, you can only integrate up to this t-value.

Answer 2

[Sam Chak]

OP: My problem is that I don't undestand how to choose the right parameters every time.

Hi @Athanasios Paraskevopoulos

In my previous post, which I have since deleted, I suggested injecting a control variable to drive the Gierer-Meinhardt activator-inhibitor system towards desired states. This approach is highly effective as it allows us to bypass concerns about the system's parameters, similar to ignoring the enemy's defense stat when launching an attack. However, upon realizing that this is an exercise provided by your professor to study the effect of parameters, I chose to remove it.

To tackle the issue of parameter selection, it is necessary to derive formulas for determining the equilibrium point. By setting both du/dt and dv/dt to zero and solving for u and v, we obtain the following formulas:

,

Ensuring positive concentrations of u and v requires satisfying the condition . This condition proves particularly helpful in your study, as it allows you to fix any three parameters {P, b, a, k} and freely choose the remaining parameter. With the singularity occurring when , this condition empowers you to select a value for the free parameter that leads to v asymptotically converging to the positive equilibrium within the vicinity of the fixed initial value of .

I believe that this represents the true essence of the study your professor wants you to grasp. It is not about blindly selecting values and creating random plots.

Example:

In the given example, {P, b, k} are fixed, and the parameter a can be freely chosen. To ensure the positivity condition is met, a must be less than 0.6750. For this particular case, the value of the free parameter a is selected as . Even if you run the simulation up to time sec, you will find that the concentration of the inhibitor v does not decrease to zero.

syms u v P b mu k a
assume(P,  {'real', 'positive'})
assume(b,  {'real', 'positive'})
assume(mu, {'real', 'positive'})
assume(k,  {'real', 'positive'})
assume(a,  {'real', 'positive'})

%% formulas to find equilibrium point
eq1     = P - a*u^2/v - mu*u == 0;
eq2     =     b*u^2   -  k*v == 0;
eqn     = [eq1, eq2];
sol     = solve(eqn, [u, v]);
usol    = simplify(sol.u, 'steps', 100)

usol = 

vsol    = simplify(sol.v, 'steps', 100)

vsol = 

%% formula to determine the parameters such that equilibrium point is positive
asol    = solve(P*b - a*k == 0, a)

asol = 

% bsol    = solve(P*b - a*k == 0, b)
% ksol    = solve(P*b - a*k == 0, k)
% Psol    = solve(P*b - a*k == 0, P)

%% parameters
value.P = 0.75;     % Production rate of u
value.b = 0.9;      % Production rate of v stimulated by u
value.mu= 0.5;      % Decay rate of u
value.k = 1;        % Decay rate of v (constant in this scenario)
value.a = 0.67;     % Ensure that a < acrit
param   = struct2array(value);

%% find equilibrium point
acrit   = double(subs(asol, [P b mu k a], param))

acrit = 0.6750

ucrit   = double(subs(usol, [P b mu k a], param))

ucrit = 0.0111

vcrit   = double(subs(vsol, [P b mu k a], param))

vcrit = 1.1111e-04

%% call ode45 solver and plot result
tspan   = [0 100];
x0      = [1.0; 0.5];   % Professor said cannot change given initial values
options = odeset('RelTol', 1e-4, 'AbsTol', 1e-6);
[t, x]  = ode45(@(t, x) ode(t, x, param), tspan, x0, options);
plot(t, x), grid on, legend('activator, u', 'inhibitor, v'), xlabel('t')
title('Responses in Gierer-Meinhardt activator-inhibitor')

%% Check if the states (u & v) settle near the critical point (equilibrium)
x(end,:)

ans = 1x2

0.0373 0.0013

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

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

%% Gierer-Meinhardt activator-inhibitor model
function dx = ode(t, x, param)
    P       = param(1); 
    b       = param(2); 
    mu      = param(3);
    k       = param(4);
    a       = param(5);
    
    u       = x(1);
    v       = x(2);
    
    dx(1,1) = P - a*u^2/v - mu*u;
    dx(2,1) =     b*u^2    - k*v;
end

[Athanasios Paraskevopoulos]

@Sam Chak Thank you very much! For your explanation and your time. I updated my code like this

%% Symbolic calculations for equilibrium points and conditions
syms u v P b mu k a
assume(P,  {'real', 'positive'})
assume(b,  {'real', 'positive'})
assume(mu, {'real', 'positive'})
assume(k,  {'real', 'positive'})
assume(a,  {'real', 'positive'})

% Equations to find equilibrium point
eq1 = P - a*u^2/v - mu*u == 0;
eq2 = b*u^2 - k*v == 0;
eqn = [eq1, eq2];
sol = solve(eqn, [u, v]);
usol = simplify(sol.u, 'steps', 100);
vsol = simplify(sol.v, 'steps', 100);

% Formula to determine the parameters such that equilibrium point is positive
asol = solve(P*b - a*k == 0, a);

%% Parameters
value.P = 0.75;     % Production rate of u
value.b = 0.9;      % Production rate of v stimulated by u
value.mu = 0.5;     % Decay rate of u
value.k = 1;        % Decay rate of v (constant in this scenario)
as_values = [0.1, 0.3,0.67]; % Different values of a to test

%% Initialize figure
figure(1);
hold on;

% Loop through each value of a and simulate using Euler method
for a_value = as_values
    value.a = a_value; % Update current value of a
    param = struct2array(value);

    %% Find equilibrium point for the current a value
    acrit = double(subs(asol, [P b mu k a], param))
    ucrit = double(subs(usol, [P b mu k a], param));
    vcrit = double(subs(vsol, [P b mu k a], param));

    %% Euler method simulation setup
    tspan = [0 100]; % Simulation time
    dt = 0.01; % Time step
    N = ceil((tspan(2) - tspan(1)) / dt); % Number of steps
    t = linspace(tspan(1), tspan(2), N);
    x = zeros(N, 2); % Initialize state variables u and v
    x(1,:) = [1.0, 0.5]; % Initial conditions

    % Euler integration loop
    for i = 1:N-1
        dx1 = value.P - value.a * x(i,1)^2 / x(i,2) - value.mu * x(i,1);
        dx2 = value.b * x(i,1)^2 - value.k * x(i,2);
        x(i+1,1) = x(i,1) + dt * dx1;
        x(i+1,2) = x(i,2) + dt * dx2;
    end

    %% Plot the results for the current a value
    plot(t, x(:,1), 'DisplayName', sprintf('u, a = %.2f', a_value));
    plot(t, x(:,2), '--', 'DisplayName', sprintf('v, a = %.2f', a_value));
end

acrit = 0.6750

acrit = 0.6750

acrit = 0.6750

%% Finalize figure with labels and legend
grid on;
legend show;
xlabel('Time t');
ylabel('Concentration');
title(' Dynamics of u and v over Time for Differen for different a values');

[Sam Chak]

Hi @Athanasios Paraskevopoulos,

I will move the computation of 'acrit' outside of the loop since we don't need to calculate the same critical value three times. The 'acrit' value will help us determine the appropriate value for 'as_values'.

%% Symbolic calculations for equilibrium points and conditions
syms u v P b mu k a
assume(P,  {'real', 'positive'})
assume(b,  {'real', 'positive'})
assume(mu, {'real', 'positive'})
assume(k,  {'real', 'positive'})
assume(a,  {'real', 'positive'})

% Equations to find equilibrium point
eq1     = P - a*u^2/v - mu*u == 0;
eq2     = b*u^2 - k*v == 0;
eqn     = [eq1, eq2];
sol     = solve(eqn, [u, v]);
usol    = simplify(sol.u, 'steps', 100);
vsol    = simplify(sol.v, 'steps', 100);

% Formula to determine the parameters such that equilibrium point is positive
asol    = solve(P*b - a*k == 0, a);

%% Parameters
value.P = 0.75;     % Production rate of u
value.b = 0.9;      % Production rate of v stimulated by u
value.mu= 0.5;      % Decay rate of u
value.k = 1;        % Decay rate of v (constant in this scenario)
param   = struct2array(value);
acrit   = double(subs(asol, [P b mu k], param))

acrit = 0.6750

%% Select values such that a < acrit
as_values = [2/9, 4/9, 6/9];    % Uniform spacing (0.2222, 0.4444, 0.6666)

%% Initialize figure
figure(1);
hold on;

% Loop through each value of a and simulate using Euler method
for j = 1:length(as_values)
    value.a = as_values(j); % Update current value of a
    param   = struct2array(value);

    %% Find equilibrium point for the current a value
    ucrit   = double(subs(usol, [P b mu k a], param))
    vcrit   = double(subs(vsol, [P b mu k a], param))

    %% Euler method simulation setup
    tspan   = [0 100];      % Simulation time
    dt      = 0.01;         % Time step
    N       = ceil((tspan(2) - tspan(1)) / dt); % Number of steps
    t       = linspace(tspan(1), tspan(2), N);
    x       = zeros(N, 2);  % Initialize state variables u and v
    x(1,:)  = [1.0, 0.5];   % Initial conditions

    % Euler integration loop
    for i = 1:N-1
        dx1     = value.P - value.a * x(i,1)^2 / x(i,2) - value.mu * x(i,1);
        dx2     = value.b * x(i,1)^2 - value.k * x(i,2);
        x(i+1,1)= x(i,1) + dt * dx1;
        x(i+1,2)= x(i,2) + dt * dx2;
    end

    %% Plot the results for the current a value
    plot(t, x(:,1),       'DisplayName', sprintf('u, a = %.4f', as_values(j)));
    plot(t, x(:,2), '--', 'DisplayName', sprintf('v, a = %.4f', as_values(j)));
end

ucrit = 1.0062

vcrit = 0.9111

ucrit = 0.5123

vcrit = 0.2362

ucrit = 0.0185

vcrit = 3.0864e-04

%% Finalize figure with labels and legend
grid on; ylim([-0.2 1.2])
legend('show', 'location', 'best');
xlabel('Time t');
ylabel('Concentration');
title('Time responses of u & v for different "a" values');