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{:});
