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

9 次查看（过去 30 天）
Athanasios Paraskevopoulos
评论： Sam Chak 2024-5-11
Ran in:
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;

请先登录，再回答此问题。

采纳的回答

Sam Chak
Sam Chak 2024-5-9
Ran in:
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
  19 个评论
显示 17更早的评论
Torsten
Torsten 2024-5-10
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.

请先登录，再进行评论。

更多回答（1 个）

Sam Chak
Sam Chak 2024-5-11
Ran in:
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
  2 个评论
Athanasios Paraskevopoulos
Athanasios Paraskevopoulos 2024-5-11
Ran in:
@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
Sam Chak 2024-5-11
Ran in:
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');

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Scenario Generation and Visualization 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by