Fixed point(equilibrium points)

How can I write a script that uses a numerical method to identify any equilibrium values(fixed points) for a system of three ordinary differential equation for a choosen set of parameters.

function Y_m3a=m3a(t,X)

%#ok<*INUSL>

global lwd n a_cca1 b_cca a_cca2 a_prr1 b_prr a_prr2 g_cca g_prr b_toc ...

g_toc k_cp k_ct k_tc k_cl k_pl k_pc k_pt

Y=zeros(1,3);

Y(1)=(((a_cca1*b_cca)+(a_cca2*b_cca*((lwd^n)/((k_cl^n) + (lwd^n)))))...

*((k_cp^n)/((k_cp^n) + (X(2)^n)))*((k_ct^n)/((k_ct^n) + (X(3)^n))))...

- (g_cca*X(1));

Y(2)=(((a_prr1*b_prr)+(a_prr2*b_prr*((lwd^n)/((k_pl^n) + (lwd^n)))))...

*((k_pc^n)/((k_pc^n) + (X(3)^n)))*((k_pt^n)/((k_pt^n) + (X(1)^n))))...

- (g_prr*X(2));

Y(3)=(b_toc*((k_tc^n)/((k_tc^n) + (X(1)^n)))) - (g_toc*X(3));

Y_m3a=Y';

end

Implementation of function

%Implementation of the m3 model from Joanito et al. 2018

% Chosen parameter see pdf note for reasons.

global lwd n a_cca1 b_cca a_cca2 a_prr1 b_prr a_prr2 g_cca g_prr b_toc ...

g_toc k_cp k_ct k_tc k_cl k_pl k_pc k_pt

lwd = 0.08324;

n = 5 ; % Hill function coefficient.

a_cca1 = 0.2042; % Fraction of cca production rate

a_cca2 = 1-a_cca1; % fraction of cca production rate associated with lwd

a_prr1 = 0.6623; % fraction of prr production rate.

a_prr2 = 1-a_prr1; % fraction of prr production rate associated with lwd

g_cca = 5.5491; % degradation rate of cca.

g_prr = 0.0105; % degradation rate of prr.

g_toc = 1.7349; % degradation rate of toc.

b_cca = g_cca; %production rate of cca

b_prr = g_prr; %production rate of prr

b_toc = g_toc; % production rate of toc

%activators and repressors coefficients used in the Hill function

k_cl = 0.1886; % lwd activating cca genes

k_cp = 0.1261; % prr genes repressing the cca genes.

k_ct = 0.1617; % toc genes repressing the cca genes.

k_pc = 0.0661; % cca genes repressing the prr genes.

k_pl = 0.0264; % lwd activating the prr genes.

k_pt = 0.1327; % toc genes repressing the prr genes.

k_tc = 0.0867; % cca genes repressing the toc genes.

y0=[0.1, 0.1, 0.1];% holds the inital condition

tspan=[0:0.01:20]; % time range with a jump of 0.01

[t,y]=ode45(@m3a,tspan,y0);% odesolver to evaluate the function and return

%t and a matrix y . The vector t will contain the

% time range

plot(t,y)

axis([0 20 0 1])

xlabel('Time')

ylabel('Concentration')

legend({'CCA1/LHY','PRR9/7','PRR5/TOC1'},'Location','east')

fun = @m3a;

x0 = [0,0,0];

Equilibrium_values = fsolve(fun,x0);

error

Not enough input arguments.

Error in m3a (line 9)

*((k_cp^n)/((k_cp^n) + (X(2)^n)))*((k_ct^n)/((k_ct^n) + (X(3)^n))))...

Error in fsolve (line 242)

fuser = feval(funfcn{3},x,varargin{:});

Answer 1

[John D'Errico]

Why would you use an ODE solver to determine equilibrium points for an ODE? Yes, I know it somehow seems logical. But this is not what ODE solvers are designed to do.

An equilibrium point is a point where the function does not change. Just set the derivatives to zero, then solve. That is, an equilibrium point is a point where Y' = 0. The same applies where you have a system of equations. As such, you use fsolve or solve or vpasolve to find that point or points, NOT an ODE solver.

[John D'Errico]

Please don't add an answer just to make a comment. Moved to a comment:

"Thank you, I just fund out that my m-file was in another folder, thats why it could not implement."