// Display mode
mode(0);

// Display warning for floating point exception
ieee(1);

// IMPLEMENTING A SIMPLE AUTO-REPRESSOR GENE NETWORK
//
// FEBRUARY 2023   M. CHAVES
//============================================================


//DELETE ALL VARIABLES FROM MEMORY
clear

//CLOSE ALL WINDOWS
xdel(winsid())


//===========================================================
// Implementation of the differential equation
function v = auto_repressor_rhs(t,x,par)

//Variables
  M = x(1);   //mRNA
  P = x(2);   //protein

//Parameters
alpha = par(1);
beta = par(2);
kP = par(3);
gM = par(4);
gP = par(5);
n = par(6);

//Initiallizing vector field
v = [];  

//Equation for mRNA
v(1) = beta * (kP^n)/(kP^n + P^n) - gM * M;
  
//Equation for protein
v(2) = alpha*M - gP * P;   

endfunction

//===========================================================
//FIND INTERSECTION VALUE FOR THE PROTEIN
function Pstst = Pintersection(P0,par)

//Parameters
alpha = par(1);
beta = par(2);
kP = par(3);
gM = par(4);
gP = par(5);
n = par(6);

Pstst = (beta/gM) * (kP^n)/(kP^n + P0^n) - (gP/alpha)*P0; 

endfunction

//========================================================
//MAIN CODE STARTS HERE

//Initial conditions, initial time, time increment, and time vector
//Order of variables: M,P

//Define parameters 
alpha = 0.2;
beta = 3.1;
kP = 10;
gM = 0.1;
gP = 0.01;
n = 2;

//Keep the same order when calling "par" in a function
par = [alpha,beta,kP,gM,gP,n];


// SOLVE SYSTEM WITH A BUILT-IN SCHEME 

x0=[10;1];
t0=0;
dt=0.05;
tvec=[t0:dt:100];

sol = ode("rk",x0,t0,tvec, list(auto_repressor_rhs,par) );

figure(1); 
plot(tvec,sol(1,:),"b-", tvec,sol(2,:),"r-");
xlabel("Time (minutes)");
ylabel("mRNA and protein");
legend("mRNA","protein");



//NUMERICALLY FIND THE STEADY STATE AS INTERSECTION BETWEEN THE TWO FUNCTIONS
Pst = fsolve(10,list(Pintersection,par));

Mst = (gP/alpha)*Pst;

ilast=length(tvec);
plot(tvec(ilast),Mst,"b*", tvec(ilast),Pst,"r*");




//PLOT THE NULLCLINES -- PHASE SPACE
Pvec = [0:0.1:3*Pst];

M1 = (beta/gM) * (kP^n)./(kP^n + Pvec.^n);
M2 = (gP/alpha)*Pvec;

figure(2);
plot(Pvec,M1,"b-", Pvec,M2,"r-");
xlabel("Protein");
ylabel("mRNA");
legend("dM/dt nullcline","dP/dt nullcline")


//PLOT ALSO THE EQUILIBRIUM POINT AND THE SOLUTION
figure(3);
plot(Pvec,M1,"b-", Pvec,M2,"r-",  Pst,Mst,'k*', sol(2,:),sol(1,:),'k--');
legend("dM/dt nullcline","dP/dt nullcline","Equilibrium point","Trajectory of solution")