function dydt = gen_net_octave(y,t)

%variable names

%proteins
%y(1) = A
%y(2) = B
%y(3) = C
%y(4) = D
%y(5) = E
%y(6) = F
%y(7) = G
%y(8) = H
%y(9) = J
%y(10) = K

%mRNAS 
%y(11) = MA
%y(12) = MB
%y(13) = MC
%y(14) = MD  
%y(15) = ME
%y(16) = MF
%y(17) = MG  
%y(18) = MH
%y(19) = MJ
%y(20) = MK

%promoters
%y(21) = PAA 
%y(22) = PBA
%y(23) = PBF 
%y(24) = PCA
%y(25) = PCD
%y(26) = PDC 
%y(27) = PDF
%y(28) = PED
%y(29) = PFEQ
%y(30) = PGC
%y(31) = PHG
%y(32) = PJK 
%y(33) = PKC

%dimers and heterodimers
%y(34) = CC
%y(35) = DD
%y(36) = FF
%y(37) = GG
%y(38) = KK
%y(39) = EQ
%y(40) = AB

%ligand
%y(41) = Q

%bound promoters
%y(42) = A-PAA
%y(43) = A-PBA
%y(44) = A-PCA
%y(45) = CC-PDC
%y(46) = CC-PGC
%y(47) = CC-PKC
%y(48) = DD-PCD
%y(49) = DD-PED
%y(50) = FF-PBF
%y(51) = FF-PDF
%y(52) = GG-PHG
%y(53) = KK-PJK
%y(54) = EQ-PFEQ

% time
%y(55) = time


%parameters for injection of ligand
%units are minutes, although the time scale is somewhat arbitrary (see readme.txt)
tstart = 1000;
tfin = 1010;
%units are molecules/cell/minute
slope = 100000;

%reaction rates and parameter values 
%units on reaction rates are molecules/cell/minute,
%although the time scale is somewhat arbitrary

%protein degradation reactions
kDA = 1;
r1= kDA*y(1);

kDB = 0.12;
r2 = kDB*y(2);

kDC = 1;
r3 = kDC*y(3);

kDegD = 1;
r4 = kDegD*y(4);

kDE = 1;
r5 = kDE*y(5);

kDF = 1;
r6 = kDF*y(6);

kDG = 1;
r7 = kDG*y(7);

kDH = 1;
r8 = kDH*y(8);

kDJ = 1;
r9 = kDJ*y(9);

kDK = 1;
r10 = kDK*y(10);

%mRNA degradation reactions
kDMA = 10;
r11 = kDMA*y(11);

kDMB = 10;
r12 = kDMB*y(12);

kDMC = 5;
r13 = kDMC*y(13);

kDMD = 5;
r14 = kDMD*y(14);

kDME = 5;
r15 = kDME*y(15);

kDMF = 5;
r16 = kDMF*y(16);

kDMG = 5;
r17 = kDMG*y(17);

kDMH = 5;
r18 = kDMH*y(18);

kDMJ = 5; 
r19 = kDMJ*y(19);

kDMK = 5;   
r20 = kDMK*y(20);

%translation reactions
kTA = 20;
r21 = kTA*y(11);

kTB = 2;
r22 = kTB*y(12);

kTC = 5;
r23 = kTC*y(13);

kTD = 5;
r24 = kTD*y(14);

kTE = 5;
r25 = kTE*y(15);

kTF = 5;
r26 = kTF*y(16);

kTG = 5;
r27 = kTG*y(17);

kTH = 5;
r28 = kTH*y(18);

kTJ = 5;
r29 = kTJ*y(19);

kTK = 5;
r30 = kTK*y(20);

%dimerization and heterodimerization reactions
kCC = 2; 
r31 = kCC*y(3)*y(3);

kDD = 2;
r32 = kDD*y(4)*y(4);

kFF = 2;
r33 = kFF*y(6)*y(6);

kGG = 2;
r34 = kGG*y(7)*y(7);

kKK = 2;
r35 = kKK*y(10)*y(10);

kEQ = 2;
r36 = kEQ*y(5)*y(41);

kAB = 2;
r37 = kAB*y(1)*y(2);


%undimerization and unheterodimerization reactions
kUC = 1;
r38 = kUC*y(34);

kUD = 1;
r39 = kUD*y(35);

kUF = 1;
r40 = kUF*y(36);

kUG = 1;
r41 = kUG*y(37);

kUK = 1;
r42 = kUK*y(38);

kUEQ = 1; 
r43 = kUEQ*y(39);

%r44 = kUAB*AB, B + 1, AB - 1 (note, as part of B&L model, AB -> B, A is degraded)
kUAB = 1;
r44 = kUAB*y(40);


%promoter binding reactions                              
kPAA = 1;
r45 = kPAA*y(1)*y(21);

kPBA = 1;
r46 = kPBA*y(1)*y(22);

kPCA = 1;
r47 = kPCA*y(1)*y(24);

kPDCC = 1;
r48 = kPDCC*y(34)*y(26);

kPGCC = 1;
r49 = kPGCC*y(34)*y(30);

kPKCC = 1;
r50 = kPKCC*y(34)*y(33);

kPCDD = 1;
r51 = kPCDD*y(35)*y(25);

kPEDD = 1;
r52 = kPEDD*y(35)*y(28);

kPBFF = 1;
r53 = kPBFF*y(36)*y(23);

kPDFF = 1;
r54 = kPDFF*y(36)*y(27);

kPHGG = 1;
r55 = kPHGG*y(37)*y(31);

kPJKK = 1;
r56 = kPJKK*y(38)*y(32);

kPFEQ = 1;
r57 = kPFEQ*y(39)*y(29);


%promoter unbinding reactions
kUPAA = 50;
r58 = kUPAA*y(42);

kUPBA = 100;
r59 = kUPBA*y(43);

kUPCA = 75;
r60 = kUPCA*y(44);

kUPDCC = 75;
r61 = kUPDCC*y(45);

kUPGCC = 75;
r62 = kUPGCC*y(46);

kUPKCC = 75;
r63 = kUPKCC*y(47);

kUPCDD = 75;
r64 = kUPCDD*y(48);

kUPEDD = 75;
r65 = kUPEDD*y(49);

kUPBFF = 75;
r66 = kUPBFF*y(50);

kUPDFF = 75;
r67 = kUPDFF*y(51);

kUPHGG = 75;
r68 = kUPHGG*y(52);

kUPJKK = 75;
r69 = kUPJKK*y(53);

kUPFEQ = 75;
r70 = kUPFEQ*y(54);

%transcription reactions
kRPAA = 50;
r71 = kRPAA*y(21);

kRAPAA = 500;
r72 = kRAPAA*y(42);

kRPBA = 0.01;
r73 = kRPBA*y(22);

kRAPBA = 50;
r74 = kRAPBA*y(43);

kRPBF = 0.005;
r75 = kRPBF*y(23);

kRFFPBF = 50;
r76 = kRFFPBF*y(50);

kRPCA = 0.01;
r77 = kRPCA*y(24);

kRAPCA =50;
r78 = kRAPCA*y(44);

kRPCD = 50;
r79 = kRPCD*y(25);

kRDDPCD = 0.01;
r80 = kRDDPCD*y(48);

kRPDC = 50;
r81 = kRPDC*y(26);

kRCCPDC = 0.01;
r82 = kRCCPDC*y(45);

kRPDF = 0.01;
r83 = kRPDF*y(27);

kRFFPDF = 50;
r84 = kRFFPDF*y(51);

kRPED = 100;
r85 = kRPED*y(28);

kRDDPED = 0.01;
r86 = kRDDPED*y(49);

kRPFEQ = 0.01;
r87 = kRPFEQ*y(29);

kREQPFEQ = 100;
r88 = kREQPFEQ*y(54);  

kRPGC = 50;
r89 = kRPGC*y(30);   

kRCCPGC = 0.01;
r90 = kRCCPGC*y(46);  

kRPHG = 50;
r91 = kRPHG*y(31);   

kRGGPHG = 0.01;
r92 = kRGGPHG*y(52);  

kRPJK = 0.01;
r93 = kRPJK*y(32);   

kRKKPJK = 50;
r94 = kRKKPJK*y(53);   

kRPKC = 50;
r95 = kRPKC*y(33);   

kRCCPKC = 0.01;
r96 = kRCCPKC*y(47);

%time 
r97 = 1;

%ligand injection 
if y(55) <= tstart
	r98 = 0;
elseif y(55) >= tfin
	r98 = 0;
else
	r98 = slope;
end;

%degradation of ligand, r99 = kDQ*Q
kDQ = 0.1;
r99 = kDQ*y(41);


%stoichiometry and ODE
%y(1) = A
dydt = [-1*r1 + r21 - r37 - r45 - r46 - r47 + r58 + r59 + r60;
%y(2) = B   
-1*r2 + r22 - r37 + r44;
%y(3) = C   
-1*r3 + r23 - 2*r31 + 2*r38;
%y(4) = D
-1*r4 + r24 - 2*r32 + 2*r39;
%y(5) = E   
-1*r5 + r25 - r36 + r43;
%y(6) = F   
-1*r6 + r26 - 2*r33 + 2*r40;
%y(7) = G   
-1*r7 + r27 - 2*r34 + 2*r41;
%y(8) = H   
-1*r8 + r28;
%y(9) = J
-1*r9 + r29;
%y(10) = K
-1*r10 + r30 - 2*r35 + 2*r42;
%y(11) = MA
-1*r11 + r71 + r72;
%y(12) = MB
-1*r12 + r73 + r74 + r75 + r76;
%y(13) = MC
-1*r13 + r77 + r78 + r79 + r80;
%y(14) = MD
-1*r14 + r81 + r82 + r83 + r84;
%y(15) = ME
-1*r15 + r85 + r86;
%y(16) = MF
-1*r16 + r87 + r88;
%y(17) = MG
-1*r17 + r89 + r90;
%y(18) = MH
-1*r18 + r91 + r92;
%y(19) = MJ
-1*r19 + r93 + r94;
%y(20) = MK
-1*r20 + r95 + r96;
%y(21) = PAA  
-1*r45 + r58;
%y(22) = PBA   
-1*r46 + r59;
%y(23) = PBF      
-1*r53 + r66;
%y(24) = PCA   
-1*r47 + r60;
%y(25) = PCD   
-1*r51 + r64;
%y(26) = PDC   
-1*r48 + r61;
%y(27) = PDF   
-1*r54 + r67;
%y(28) = PED   
-1*r52 + r65;
%y(29) = PFEQ  
-1*r57 + r70;
%y(30) = PGC   
-1*r49 + r62;
%y(31) = PHG
-1*r55 + r68;
%y(32) = PJK
-1*r56 + r69;
%y(33) = PKC
-1*r50 + r63;
%y(34) = CC
r31 - r38 - r48 - r49 - r50 + r61 + r62 + r63;
%y(35) = DD
r32 - r39 - r51 - r52 + r64 + r65;
%y(36) = FF 
r33 - r40 - r53 - r54 + r66 + r67;
%y(37) = GG
r34 - r41 - r55 + r68;
%y(38) = KK
r35 - r42 - r56 + r69;
%y(39) = EQ
r36 - r43 - r57 + r70;
%y(40) = AB
r37 - r44;
%y(41) = Q 
-r36 + r43 + r98 - r99;
%y(42) = A-PAA
r45 - r58;
%y(43) = A-PBA
r46 - r59;
%y(44) = A-PCA
r47 - r60;
%y(45) = CC-PDC
r48 - r61;
%y(46) = CC-PGC
r49 - r62;
%y(47) = CC-PKC
r50 - r63;
%y(48) = DD-PCD
r51 - r64;
%y(49) = DD-PED
r52 - r65;
%y(50) = FF-PBF
r53 - r66;
%y(51) = FF-PDF
r54 - r67;
%y(52) = GG-PHG
r55 - r68;
%y(53) = KK-PJK
r56 - r69;
%y(54) = EQ-PFEQ
r57 - r70;
%y(55) = time
r97];

endfunction