function x_output = E2depICIpalbo_model(~,x,PAR,ICI,palbo,abema)
% variables
x(x < 0) = 0;
E2media = x(1);
E2cell = x(2);
ER = x(3);
E2ER = x(4);
E2NSB = x(5);
ICIER = x(6);
rescyclinD1palbo = x(7);
cyclinD1 = x(8);
cdk46 = x(9);
cdk46palbo = x(10);
cdk46abema = x(11);
cyclinD1cdk46 = x(12);
cyclinD1cdk46p21 = x(13);
cyclinD1cdk46palbo = x(14);
cyclinD1cdk46abema = x(15);
cyclinD1cdk46p21palbo = x(16);
cyclinD1cdk46p21abema = x(17);
cMyc = x(18);
p21 = x(19);
cyclinE = x(20);
cyclinEp21 = x(21);
Rb = x(22);
pRb = x(23);
ppRb = x(24);
respropalbo = x(25);
Nalive = x(26);
Ndead = x(27);

% Parameter
% (1).Diffusion rate of E2 
k_diff = PAR(1);

% (2).Binding rate between non-specific binding and E2
kb_NSB = PAR(2);

% (3).Unbinding rate between non-specific binding and E2
kub_NSB = PAR(3);

% (4).Volume of MCF7 cell
Vol_1cell = PAR(4);

% (5).Volume of media
Vol_media = PAR(5);

% (6).Translation rate of ER
k_ER = PAR(6);

% (7).Degradation rate of ER
kd_ER = PAR(7);

% (8).Degradation rate of E2ER
kd_E2ER = PAR(8);

% (9).Binding rate between E2cell and ER
kb_E2ER = PAR(9);

% (10).Unbinding rate between E2cell and ER
kub_E2ER = PAR(10);

% (11).Binding rate between ICI and ER
kb_ICIER = PAR(11);

% (12).Unbinding rate between ICI and ER
kub_ICIER = PAR(12);

% (13).Degradation rate of ICIER
kd_ICIER = PAR(13);

% (14).Translation rate of cyclinD1
k_cyclinD1 = PAR(14);

% (15).Degradation rate of cyclinD1
kd_cyclinD1 = PAR(15);

% (16).Increased cyclinD1 translation by E2ER
k_cyclinD1E2ER = PAR(16);

% (17).Parameter 1 of cyclinD1 increased translation by E2ER
p_cyclinD1E2ER_1 = PAR(17);

% (18).Parameter 2 of cyclinD1 increased translation by E2ER
p_cyclinD1E2ER_2 = PAR(18);

% (19).Binding rate between cyclinD1 and cdk46
kb_cyclinD1cdk46 = PAR(19);

% (20).Unbinding rate between cyclinD1 and cdk46
kub_cyclinD1cdk46 = PAR(20);

% (21).Generation rate of rescyclinD1palbo
k_rescyclinD1palbo = PAR(21);

% (22).Parameter 1 of rescyclinD1palbo increased by palbo
p_rescyclinD1palbo_1 = PAR(22);

% (23).Parameter 2 of rescyclinD1palbo increased by palbo
p_rescyclinD1palbo_2 = PAR(23);

% (24).Degradation rate of rescyclinD1palbo
kd_rescyclinD1palbo = PAR(24);

% (25).Increase rate of cyclinD1 by rescyclinD1palbo
k_cyclinD1palbo = PAR(25);

% (26).Parameter 1 of cyclinD1 increased by rescyclinD1palbo
p_cyclinD1palbo_1 = PAR(26);

% (27).Parameter 2 of cyclinD1 increased by rescyclinD1palbo
p_cyclinD1palbo_2 = PAR(27);

% (28).Translation rate of cdk46
k_cdk46 = PAR(28);

% (29).Degradation rate of cdk46
kd_cdk46 = PAR(29);

% (30).Binding rate between cdk46 and palbo
kb_cdk46palbo = PAR(30);

% (31).Unbinding rate between cdk46 and palbo
kub_cdk46palbo = PAR(31);

% (32) Binding rate between cdk46 and abema
kb_cdk46abema = PAR(32);

% (33).Unbinding rate between cdk46 and abema
kub_cdk46abema = PAR(33);

% (34).Degradation rate of cyclinD1cdk46
kd_cyclinD1cdk46 = PAR(34);

% (35).Binding rate between cyclinD1cdk46 and palbo
kb_cyclinD1cdk46palbo = PAR(35);

% (36).Unbinding rate between cyclinD1cdk46 and palbo
kub_cyclinD1cdk46palbo = PAR(36);

% (37).Binding rate between cyclinD1cdk46 and abema
kb_cyclinD1cdk46abema = PAR(37);

% (38).Unbinding rate between cyclinD1cdk46 and abema
kub_cyclinD1cdk46abema = PAR(38);

% (39).Bindinding rate beteween cyclinD1cdk46 and p21
kb_cyclinD1cdk46p21 = PAR(39);

% (40).Unbinding rate between cyclinD1cdk46 and p21
kub_cyclinD1cdk46p21 = PAR(40);

% (41).Degradation rate of cyclinD1cdk46p21
kd_cyclinD1cdk46p21 = PAR(41);

% (42).Binding rate between cyclinD1cdk46p21 and palbo
kb_cyclinD1cdk46p21palbo = PAR(42);

% (43).Unbinding rate between cyclinD1cdk46p21 and palbo
kub_cyclinD1cdk46p21palbo = PAR(43); 

% (44).Binding rate between cyclinD1cdk46p21 and abema
kb_cyclinD1cdk46p21abema = PAR(44);

% (45).Unbinding rate between cyclinD1cdk46p21 and abema
kub_cyclinD1cdk46p21abema = PAR(45);

% (46).Binding rate between cyclinD1cdk46palbo and p21
kb_cyclinD1cdk46palbop21 = PAR(46);

% (47).Unbinding rate between cyclinD1cdk46palbo and p21
kub_cyclinD1cdk46palbop21 = PAR(47);

% (48).Binding rate between cyclinD1cdk46abema and p21
kb_cyclinD1cdk46abemap21 = PAR(48);

% (49).Unbindng rate between cyclinD1cdk46abema and p21
kub_cyclinD1cdk46abemap21 = PAR(49);

% (50).Degradation rate of cyclinD1cdk46palbo or cyclinD1cdk46abema
kd_cyclinD1cdk46palboabema = PAR(50);

% (51).Translation rate of cMyc
k_cMyc = PAR(51);

% (52).Degradation rate of cMyc
kd_cMyc = PAR(52);

% (53).Increased translation of cMyc by E2ER
k_cMycE2ER = PAR(53);

% (54).Parameter 1 of cMyc increased translation by E2ER
p_cMycE2ER_1 = PAR(54);

% (55).Parameter 2 of cMyc increased translation by E2ER
p_cMycE2ER_2 = PAR(55);

% (56).Increased translation of cMyc by ppRb
k_cMycppRb = PAR(56);

% (57).Parameter 1 of cMyc increased translation by ppRb
p_cMycppRb_1 = PAR(57);

% (58).Parameter 2 of cMyc increased translation by ppRb
p_cMycppRb_2 = PAR(58);

% (59).Translation rate of p21
k_p21 = PAR(59);

% (60).Degradation rate of p21
kd_p21 = PAR(60);

% (61).Rate of p21 translation inhibited by cMyc
k_p21cMyc = PAR(61);

% (62).Parameter 1 of p21 inhibited translation by cMyc
p_p21cMyc_1 = PAR(62);

% (63).Parameter 2 of p21 inhibitied translation by cMyc
p_p21cMyc_2 = PAR(63);

% (64).Translation rate of cyclinE
k_cyclinE = PAR(64);

% (65).Degradation rate of cyclinE
kd_cyclinE = PAR(65);

% (66).Increased translation rate of cyclinE by E2ER
k_cyclinEE2ER = PAR(66);

% (67).Parameter 1 of cyclinE increased translation by E2ER
p_cyclinEE2ER_1 = PAR(67);

% (68).Parameter 2 of cyclinE increased translation by E2ER
p_cyclinEE2ER_2 = PAR(68);

% (69).Binding rate between cyclinE and p21
kb_cyclinEp21 = PAR(69);

% (70).Unbinding rate between cyclinE and p21
kub_cyclinEp21 = PAR(70);

% (71).Translation rate of Rb
k_Rb = PAR(71);

% (72).Degradation rate of Rb
kd_Rb = PAR(72);

% (73).Increased Rb translation by ppRb
k_RbppRb = PAR(73);

% (74).Parameter 1 of Rb increased translation by ppRb
p_RbppRb_1 = PAR(74);

% (75).Parameter 2 of Rb increased translation by ppRb
p_RbppRb_2 = PAR(75);

% (76).Phosphorylation rate of Rb by cyclinD1cdk46
k_RbcyclinD1cdk46 = PAR(76);

% (77).Dephosphorylation rate of pRb
k_pRbdepho = PAR(77);

% (78).Degradation rate of pRb
kd_pRb = PAR(78);

% (79).Phosphorylation rate of pRb by cyclinE
k_pRbcyclinE = PAR(79);

% (80).Dephosphorylation rate of ppRb
k_ppRbdepho = PAR(80);

% (81).Degradation rate of ppRb
kd_ppRb = PAR(81);

% (82).Parameter 1 of cyclinD1cdk46 kinase activity
p_cyclinD1cdk46_1 = PAR(82);

% (83).Parameter 2 of cyclinD1cdk46 kinase activity
p_cyclinD1cdk46_2 = PAR(83);

% (84).Parameter 1 of Rb phosphorylation by cyclinD1cdk46
p_cyclinD1cdk46Rb_1 = PAR(84);

% (85).Parameter 2 of Rb phosphorylation by cyclinD1cdk46
p_cyclinD1cdk46Rb_2 = PAR(85);

% (86).Parameter 1 of pRb dephosphorylation
p_pRb_1 = PAR(86);

% (87).Parameter 2 of pRb dephosphorylation
p_pRb_2 = PAR(87);

% (88).Parameter 1 of pRb phosphorylation by cyclinE
p_cyclinEpRb_1 = PAR(88);

% (89).Parameter 2 of pRb phosphorylation by cyclinE
p_cyclinEpRb_2 = PAR(89);

% (90).Parameter 1 of ppRb dephosphorylation
p_ppRb_1 = PAR(90);

% (91).Parameter 2 of ppRb dephosphorylation
p_ppRb_2 = PAR(91);

% (92).Basal proliferation rate
k_pro = PAR(92);

% (93).Proliferation rate increased by ppRb
k_proppRb = PAR(93);

% (94).Parameter 1 of proliferation increased by ppRb
p_proppRb_1 = PAR(94);

% (95).Parameter 2 of proliferation increased by ppRb
p_proppRb_2 = PAR(95);

% (96).Generation rate of respropalbo
k_respropalbo = PAR(96);

% (97).Parameter 1 of increasing rate of respropalbo by palbo
p_respropalbo_1 = PAR(97);

% (98).Parameter 2 of increasing rate of respropalbo by palbo
p_respropalbo_2 = PAR(98);

% (99).Parameter of degradation rate of respropalbo
p_respropalbo_kd = PAR(99);

% (100).Degradation rate of respropalbo
kd_respropalbo = PAR(100);

% (101).Rate of proliferation inhibited by respropalbo
k_propalbo = PAR(101);

% (102).Parameter 1 of proliferation inhibited by respropalbo
p_propalbo_1 = PAR(102);

% (103).Parameter 2 of proliferation inhibited by respropalbo
p_propalbo_2 = PAR(103);

% (104).Carrying capacity
k_carrying = PAR(104);

% (105).Basal death rate
k_death = PAR(105);

% (106).Lysis rate of dead cell
k_lysis = PAR(106);

% Model

% Total number of cells equals number of alive cells plus number of dead cells
N = Nalive + Ndead; 

dNalive = k_pro * (1 + k_proppRb * ppRb^p_proppRb_2/(p_proppRb_1^p_proppRb_2 + ppRb^p_proppRb_2))/(1 + k_propalbo * respropalbo^p_propalbo_2/(p_propalbo_1^p_propalbo_2 + respropalbo^p_propalbo_2))...
    * Nalive * (1 - N/k_carrying)... % Basal proliferation, increased proliferation by ppRb and inhibited proliferation by respropalbo
    - k_death * Nalive; % Basal death

dNdead = k_death * Nalive... % Basal death
    - k_lysis * Ndead; % Lysis of dead cells

% Hill function for respropalbo
hillfun_respropalbo = palbo^p_respropalbo_2/(palbo^p_respropalbo_2 + p_respropalbo_1^p_respropalbo_2);

drespropalbo = k_respropalbo * hillfun_respropalbo... % Generation of respropalbo by palbo
    - kd_respropalbo*respropalbo/(1 + p_respropalbo_kd * hillfun_respropalbo); % Degradation of respropalbo 

% E2 concentration changes in media
dE2media = k_diff * (N * Vol_1cell)/Vol_media * (E2cell - E2media);

dE2cell = - k_diff * (E2cell - E2media)- (dNalive + dNdead)/N * E2cell... % E2 concentration changes in cell
    - kb_E2ER * E2cell * ER + kub_E2ER * E2ER... % Binding and unbinding between ER and E2cell
    - kb_NSB * E2cell + kub_NSB * E2NSB... % Binding and unbinding between non-specific binding and E2cell in the cell
    + kd_E2ER * E2ER; % Degradation of E2ER

dER = k_ER - kd_ER * ER... % Translation and degradation of ER
    - kb_E2ER * E2cell * ER + kub_E2ER * E2ER... % Binding and unbinding between ER and E2cell
    - kb_ICIER * ICI * ER + kub_ICIER * ICIER; % Binding and unbinding between ER and ICI

dE2ER = - kd_E2ER * E2ER... % Degradation of E2ER
    + kb_E2ER * E2cell * ER - kub_E2ER * E2ER; % Binding and unbinding between ER and E2cell

dE2NSB = kb_NSB *  E2cell - kub_NSB * E2NSB; % Binding and unbinding between non-specific binding and E2cell

dICIER = kb_ICIER * ICI * ER - kub_ICIER * ICIER... % Binding and unbinding between ICI and ER
    - kd_ICIER * ICIER; % Degradation of ICIER


hillfun_rescyclinD1palbo = palbo^p_rescyclinD1palbo_2/(palbo^p_rescyclinD1palbo_2 + p_rescyclinD1palbo_1^p_rescyclinD1palbo_2);

drescyclinD1palbo = k_rescyclinD1palbo * hillfun_rescyclinD1palbo... % Generation of rescyclinD1palbo by palbo
    - kd_rescyclinD1palbo * rescyclinD1palbo; % Degradation of rescyclinD1palbo

dcyclinD1 = - kd_cyclinD1 * cyclinD1... % Degradation of cyclinD1
    + k_cyclinD1 * (1 + k_cyclinD1E2ER * E2ER^p_cyclinD1E2ER_2/(p_cyclinD1E2ER_1^p_cyclinD1E2ER_2 + E2ER^p_cyclinD1E2ER_2))...% Basal translation of cyclinD1 and the increased translation by E2ER
    - kb_cyclinD1cdk46 * cyclinD1 * cdk46 + kub_cyclinD1cdk46 * cyclinD1cdk46... % Binding and unbinding between cyclinD1 and cdk46
    + k_cyclinD1palbo * rescyclinD1palbo^p_cyclinD1palbo_2/(p_cyclinD1palbo_1^p_cyclinD1palbo_2 + rescyclinD1palbo^p_cyclinD1palbo_2);% Increased translation of cyclinD1 by rescyclinD1palbo

dcdk46 = k_cdk46 - kd_cdk46 * cdk46... % Translation and degradation of cdk46
    - kb_cyclinD1cdk46 * cyclinD1 * cdk46 + kub_cyclinD1cdk46 * cyclinD1cdk46...% Binding and unbinding between cyclinD1 and cdk46
    - kb_cdk46palbo * cdk46 * palbo + kub_cdk46palbo * cdk46palbo...  % Binding and unbinding between palbo and cdk46
    - kb_cdk46abema * cdk46 * abema + kub_cdk46abema * cdk46abema; % Binding and unbinding between abema and cdk46
      
dcdk46palbo = - kd_cdk46 * cdk46palbo... % Degradation of cdk46palbo
    + kb_cdk46palbo * cdk46 * palbo - kub_cdk46palbo * cdk46palbo; % Binding and unbinding between palbo and cdk46

dcdk46abema = - kd_cdk46 * cdk46abema... % Degradation of cdk46abema
    + kb_cdk46abema * cdk46 * abema - kub_cdk46abema * cdk46abema; % Binding and unbinding between abema and cdk46

dcyclinD1cdk46 = - kd_cyclinD1cdk46 * cyclinD1cdk46...% Degradation of cyclinD1cdk46
     + kb_cyclinD1cdk46 * cyclinD1 * cdk46 - kub_cyclinD1cdk46 * cyclinD1cdk46... % Binding and unbinding between cyclinD1 and cdk46
     - kb_cyclinD1cdk46p21 * cyclinD1cdk46 * p21 + kub_cyclinD1cdk46p21 * cyclinD1cdk46p21... % Binding and unbinding between p21 and cyclinD1cdk46
     - kb_cyclinD1cdk46palbo * cyclinD1cdk46 * palbo + kub_cyclinD1cdk46palbo * cyclinD1cdk46palbo... % Binding between palbo and cyclinD1cdk46
     - kb_cyclinD1cdk46abema * cyclinD1cdk46 * abema + kub_cyclinD1cdk46abema * cyclinD1cdk46abema; % Binding between abema and cyclinD1cdk46
 
dcyclinD1cdk46p21 = - kd_cyclinD1cdk46p21 * cyclinD1cdk46p21... % Degradation of cyclinD1cdk46p21
    + kb_cyclinD1cdk46p21 * cyclinD1cdk46 * p21 - kub_cyclinD1cdk46p21 * cyclinD1cdk46p21... % Binding and unbinding between p21 and cyclinD1cdk46
    - kb_cyclinD1cdk46p21palbo * cyclinD1cdk46p21 * palbo + kub_cyclinD1cdk46p21palbo * cyclinD1cdk46p21palbo... % Binding and unbinding between palbo and cyclinD1cdk46p21
    - kb_cyclinD1cdk46p21abema * cyclinD1cdk46p21 * abema + kub_cyclinD1cdk46p21abema * cyclinD1cdk46p21abema; % Binding and unbinding between abema and cyclinD1cdk46p21

dcyclinD1cdk46palbo = - kd_cyclinD1cdk46palboabema * cyclinD1cdk46palbo... % Degradation of cyclinD1cdk46palbo
    + kb_cyclinD1cdk46palbo * cyclinD1cdk46 * palbo - kub_cyclinD1cdk46palbo * cyclinD1cdk46palbo... % Binding and unbinding between palbo and cyclinD1cdk46
    - kb_cyclinD1cdk46palbop21 * cyclinD1cdk46palbo * p21 + kub_cyclinD1cdk46palbop21 * cyclinD1cdk46p21palbo; % Binding and unbinding between p21 and cyclinD1cdk46palbo

dcyclinD1cdk46abema = - kd_cyclinD1cdk46palboabema * cyclinD1cdk46abema... % Degradation of cyclinD1cdk46abema
    + kb_cyclinD1cdk46abema * cyclinD1cdk46 * abema - kub_cyclinD1cdk46abema * cyclinD1cdk46abema... % Binding and unbinding between abema and cyclinD1cdk46
    - kb_cyclinD1cdk46abemap21 * cyclinD1cdk46abema * p21 + kub_cyclinD1cdk46abemap21 * cyclinD1cdk46p21abema; % Binding and unbinding between p21 and cyclinD1cdk46abema

dcyclinD1cdk46p21palbo = - kd_cyclinD1cdk46p21 * cyclinD1cdk46p21palbo... % Degradation of cyclinD1cdk46p21palbo
    + kb_cyclinD1cdk46p21palbo * cyclinD1cdk46p21 * palbo - kub_cyclinD1cdk46p21palbo * cyclinD1cdk46p21palbo... % Binding and unbinding between palbo and cyclinD1cdk46p21
    + kb_cyclinD1cdk46palbop21 * cyclinD1cdk46palbo * p21 - kub_cyclinD1cdk46palbop21 * cyclinD1cdk46p21palbo; % Binding and unbinding between p21 and cyclinD1cdk46palbo

dcyclinD1cdk46p21abema = - kd_cyclinD1cdk46p21 * cyclinD1cdk46p21abema... % Degradation of cyclinD1cdk46p21abema
    + kb_cyclinD1cdk46p21abema * cyclinD1cdk46p21 * abema - kub_cyclinD1cdk46p21abema * cyclinD1cdk46p21abema... % Binding and unbinding between abema and cyclinD1cdk46p21
    + kb_cyclinD1cdk46abemap21 * cyclinD1cdk46abema * p21 - kub_cyclinD1cdk46abemap21 * cyclinD1cdk46p21abema; % Binding and unbinding between p21 and cyclinD1cdk46abema

dcMyc = - kd_cMyc * cMyc... % Degradation of cMyc
    + k_cMyc * (1 + k_cMycE2ER * E2ER^p_cMycE2ER_2/(p_cMycE2ER_1^p_cMycE2ER_2 + E2ER^p_cMycE2ER_2)... % Basal translation of cMyc and the increased translation by E2ER
    + k_cMycppRb * ppRb^p_cMycppRb_2/(p_cMycppRb_1^p_cMycppRb_2 + ppRb^p_cMycppRb_2)); % Increased translation of cMyc by ppRb

dp21 = k_p21 - kd_p21 * p21... % % Translation and degradation of p21
    + k_p21cMyc * (p_p21cMyc_1^p_p21cMyc_2/(p_p21cMyc_1^p_p21cMyc_2 + cMyc^p_p21cMyc_2))... % Inhibition of translation by cMyc
    - kb_cyclinD1cdk46p21 * cyclinD1cdk46 * p21 + kub_cyclinD1cdk46p21 * cyclinD1cdk46p21... % Binding and unbinding between p21 and cyclinD1cdk46
    - kb_cyclinD1cdk46palbop21 * cyclinD1cdk46palbo * p21 + kub_cyclinD1cdk46palbop21 * cyclinD1cdk46p21palbo... % Binding and unbinding between p21 and cyclinD1cdk46palbo
    - kb_cyclinD1cdk46abemap21 * cyclinD1cdk46abema * p21 + kub_cyclinD1cdk46abemap21 * cyclinD1cdk46p21abema... % Binding and unbinding between p21 and cyclinD1cdk46abema
    - kb_cyclinEp21 * cyclinE * p21 + kub_cyclinEp21 * cyclinEp21; % Binding and unbinding between p21 and cyclinE

dcyclinE = - kd_cyclinE * cyclinE... % Degradation of cyclinE 
    + k_cyclinE * (1 + k_cyclinEE2ER * E2ER^p_cyclinEE2ER_2/(p_cyclinEE2ER_1^p_cyclinEE2ER_2 + E2ER^p_cyclinEE2ER_2))... % Basal translation of cyclinE and the increased translation by E2ER
    - kb_cyclinEp21 * cyclinE * p21 + kub_cyclinEp21 * cyclinEp21; % Binding and unbinding between p21 and cyclinE
    
dcyclinEp21 = - kd_cyclinE * cyclinEp21... % Degradation of cyclinEp21
    + kb_cyclinEp21 * cyclinE * p21 - kub_cyclinEp21 * cyclinEp21; % Binding and unbinding between p21 and cyclinE


dRb = k_Rb - kd_Rb * Rb...  % Degradation and basal translation of Rb
    + k_RbppRb * ppRb^p_RbppRb_2/(p_RbppRb_1^p_RbppRb_2 + ppRb^p_RbppRb_2)...% Increased translation of Rb by ppRb
    - k_RbcyclinD1cdk46 * cyclinD1cdk46^p_cyclinD1cdk46_2/(p_cyclinD1cdk46_1^p_cyclinD1cdk46_2 + cyclinD1cdk46^p_cyclinD1cdk46_2)...
    * Rb^p_cyclinD1cdk46Rb_2/(p_cyclinD1cdk46Rb_1^p_cyclinD1cdk46Rb_2 + Rb^p_cyclinD1cdk46Rb_2)... % Phosphorylation of Rb by cyclinD1cdk46
    + k_pRbdepho * pRb^p_pRb_2/(p_pRb_1^p_pRb_2 + pRb^p_pRb_2); % Dephosphorylation of pRb

dpRb = - kd_pRb * pRb... % Degradation of pRb
    + k_RbcyclinD1cdk46 * cyclinD1cdk46^p_cyclinD1cdk46_2/(p_cyclinD1cdk46_1^p_cyclinD1cdk46_2 + cyclinD1cdk46^p_cyclinD1cdk46_2)...
    * Rb^p_cyclinD1cdk46Rb_2/(p_cyclinD1cdk46Rb_1^p_cyclinD1cdk46Rb_2 + Rb^p_cyclinD1cdk46Rb_2)... % Phosphorylation of Rb by cyclinD1cdk46
    - k_pRbdepho * pRb^p_pRb_2/(p_pRb_1^p_pRb_2 + pRb^p_pRb_2)... % Dephosphorylation of pRb
    - k_pRbcyclinE * cyclinE * pRb^p_cyclinEpRb_2/(p_cyclinEpRb_1^p_cyclinEpRb_2 + pRb^p_cyclinEpRb_2)... % Phosphorylation of pRb by cyclinE
    + k_ppRbdepho * ppRb^p_ppRb_2/(p_ppRb_1^p_ppRb_2 + ppRb^p_ppRb_2); % Dephosphorylation of ppRb

dppRb = - kd_ppRb * ppRb... % Degradation of ppRb
    + k_pRbcyclinE * cyclinE * pRb^p_cyclinEpRb_2/(p_cyclinEpRb_1^p_cyclinEpRb_2 + pRb^p_cyclinEpRb_2)... % Phosphorylation of pRb by cyclinE
    - k_ppRbdepho * ppRb^p_ppRb_2/(p_ppRb_1^p_ppRb_2 + ppRb^p_ppRb_2); % Dephosphorylation of ppRb

x_output(1,:) = dE2media;
x_output(2,:) = dE2cell;
x_output(3,:) = dER;
x_output(4,:) = dE2ER;
x_output(5,:) = dE2NSB;
x_output(6,:) = dICIER;
x_output(7,:) = drescyclinD1palbo;
x_output(8,:) = dcyclinD1;
x_output(9,:) = dcdk46;
x_output(10,:) = dcdk46palbo;
x_output(11,:) = dcdk46abema;
x_output(12,:) = dcyclinD1cdk46;
x_output(13,:) = dcyclinD1cdk46p21;
x_output(14,:) = dcyclinD1cdk46palbo;
x_output(15,:) = dcyclinD1cdk46abema;
x_output(16,:) = dcyclinD1cdk46p21palbo;
x_output(17,:) = dcyclinD1cdk46p21abema;
x_output(18,:) = dcMyc;
x_output(19,:) = dp21;
x_output(20,:) = dcyclinE;
x_output(21,:) = dcyclinEp21;
x_output(22,:) = dRb;
x_output(23,:) = dpRb;
x_output(24,:) = dppRb;
x_output(25,:) = drespropalbo;
x_output(26,:) = dNalive;
x_output(27,:) = dNdead;

end