FN_ODE_PID_LinTX.m

function dYdt = FN_ODE_PID_LinTX(t,y)
    load Par_ODE.mat
	
	Ys  = k.Y;
    
    %% Perturbation
    k.(k.Pn) = k.(k.Pn)*k.P(2,sum([t>=k.P(1,:)]));
	
	Y   = k.Y - Ys;

    %% Species:
    z1   = y(1);
    z2   = y(2);
    x1   = y(3);
    xC   = y(4);
    a    = y(5);
    m    = y(6);
    I    = y(7);

    %% Steady states:
    Z1s  = k.SS(1);
    Z2s  = k.SS(2);
    X1s  = k.SS(3);
    XCs  = k.SS(4);
    As   = k.SS(5);
    Ms   = k.SS(6);
    
    %% Events/functions:
    O0th = @(c)     [c];
    O1th = @(X,c)   [c*X];
    O2th = @(X,Y,c) [c*X*Y];
    OMMp = @(X,c,d) [c*X/(X+d)];
    OMMn = @(X,c,d) [c*d/(X+d)];
    
    OMLn = @(X,c,d) [c*d/((X+d)^2)];
    
    %% ODEs
    F0 = [O0th(k.mu*Ys)  - O2th(Z1s,Z2s,k.et)    - O1th(Z1s,k.gD);
          O1th(XCs,k.th) - O2th(Z1s,Z2s,k.et)    - O1th(Z2s,k.gD);
          O1th(Z1s,k.bI)  ...
            + OMMn(k.th*XCs,k.bP*Ys,k.mu*Ys) + O1th(As,k.bD) ...
                             - O1th(X1s,k.g1)    - O1th(X1s,k.gD);
          O1th(X1s,k.bC)  - O1th(XCs,k.gC)       - O1th(XCs,k.gD);
          O1th(Ms ,k.bA)  - OMMp(As,XCs*k.gA,k.KA) ...
                             - O1th(As,k.gA0)    - O1th(As,k.gD);
          O0th(k.bM*Ys)  - OMMp(Ms,As*k.gM,k.KM) - O1th(Ms,k.gD);
          0];
      
    kI   = k.bI*Z1s/(Z1s+Z2s);
    kP   = (1/k.th)*(k.bP*Ys*k.mu*Ys/((XCs+(k.mu*Ys))^2));
    kD   = (1/k.th)*(k.bD*k.gA/(k.bA*k.gM));
    
    e    = O0th(k.mu*Y)  - O1th(xC,k.th);
    dXC  = O1th(x1,k.bC) - O1th(xC,k.gC) - O1th(xC,k.gD);
			
    dYdt = F0 + ...
	       [O0th(k.mu*Y)   - O1th(z1,k.et*Z2s) - O1th(z2,k.et*Z1s) ...
												 - O1th(z1,k.gD);
            O1th(xC,k.th)  - O1th(z1,k.et*Z2s) - O1th(z2,k.et*Z1s) ...
												 - O1th(z2,k.gD);
            O1th(I+Z1s,kI) + O1th(e,kP) - O1th(dXC,kD*k.th) ...
                              - O1th(x1,k.g1)    - O1th(x1,k.gD);
            O1th(x1,k.bC)  - O1th(xC,k.gC)       - O1th(xC,k.gD);
            O1th(m ,k.bA) ...
			  - OMMp(As,xC*k.gA,k.KA) - OMLn(As,a*XCs*k.gA,k.KA) ...
                              - O1th(a,k.gA0)    - O1th(a ,k.gD);
            O0th(k.bM*Y) ...
			  - OMMp(Ms,a*k.gM,k.KM) - OMLn(Ms,m*As*k.gM,k.KM) ...
												 - O1th(m ,k.gD);
            e];