function [dxdt] = mi_lux(x,param)

  % INDUCIBLE MODEL LuxR/LuxI: 
    % Establishes model equations for the genetic circuit LuxR/LuxI of 
    % the designed circuit. It gives us back the variation of the 
    % concentrations of the chemical species of the reactions in each 
    % instant of time. 
    
  % Abbreviations:

    % R: LuxR protein
    % I: protein of interest (such as GFP)
    % R·A: (LuxR·AHL)
    % (R·A)2: (LuxR·AHL) dimer
    % A: AHL molecule
    % mRNA: messenger RNA

    % param: parameters vector
        % pI: mGFP translation rate (min-1)
        % kI: gGFP transcription rate (molec.min-1)
        % pN_luxI: plasmid number of LuxI (plasmids)
        % dmI: mGFP degradation rate (min-1)
        % kdLux: Dissociation constant of (R·A)2 to the PluxR promoter (molec)
        % alphaI: PluxR promoter basal activity (adimensional)
        % dI: GFP degradation rate (min-1)
        % k_1: Dissociation rate of (R·A) (min-1)
        % kd1: Dissociation constant of (R·A) (molec)
        % pR: mR translation rate (min-1)
        % cR: Constitutive transcription rate of LuxR protein (min-1)
        % dmR: mR degradation rate (min-1)
        % k_2: Dissociation rate of (R·A)2 (min-1)
        % kd2: Dissociation constant of (R·A)2 (molec)
        % dRA: (R·A) degradation rate (min-1)
        % dR: R degradation rate (min-1)
        % dRA2: (R·A)2 degradation rate (min-1)
        % D: Diffusion coefficient of AHL through the cell membrane (min-1)
        % Vc: Vcell/Vext ratio (min-1)
        % dA: A degradation rate (min-1)
        % dAe: Ae degradation rate (min-1)
        % N: number of cells (cells)

n1 = x(1);
n2 = x(2);
n3 = x(3);
n4 = x(4);
n5 = x(5);

pI = param(1);
kI = param(2);
pN_luxI = param(3);
dmI = param(4);
kdLux = param(5);
alphaI = param(6);
dI = param(7);
k_1 = param(8);
kd1 = param(9);
pR = param(10);
cR = param(11);
dmR = param(12);
k_2 = param(13);
kd2 = param(14);
dRA = param(15);
dR = param(16);
dRA2 = param(17);
D = param(18);
Vc = param(19);
dA = param(20);
dAe = param(21);
N = param(22);

% n1: LuxI / Protein of interest
c1 = (pI*kI *pN_luxI)./dmI;
c3 = kdLux + n3; 
dxdt1 = c1*( kdLux/c3 + (alphaI*n3)/c3 ) - dI*n1; %% Represor
% dxdt1 = c1*( alphaI*kdLux/c3 + n3/c3 ) - dI*n1; %% Activator

% n2: LuxR
c2 = k_1./kd1;
c4 = (pR*cR)./dmR;

% Algebraic equation for (LuxR.AHL)
c6 = (2*k_2*kd1*n3 + k_1*n2*n4);
c7 = (8*k_2)*c6;
c8 = k_1 + dRA;
c9 = (kd1*kd2)*c8^2;
c10 = (kd2*c8)/(4*k_2);
c11 = c7/c9 + 1;
MO  = c10*(sqrt(abs(c11)) - 1);  % This is LuxR.AHL
assignin('base', 'M', MO )

dxdt2 = c4 + k_1*MO  - (c2*n4 + dR)*n2; 

% n3: (LuxR.AHL)2
c5 = k_2./kd2;
dxdt3 = c5*MO^2 - (k_2 + dRA2)*n3;

% n4: AHLint
dxdt4 = k_1*MO  + D*Vc*n5  - (D + c2*n2 + dA)*n4;

% n5: AHLext
dxdt5 = D*N*n4 - (D*Vc*N + dAe)*n5; 


dxdt = [dxdt1;dxdt2;dxdt3;dxdt4;dxdt5];

end