% Simple model of a genetic switch

% This function implements the basic model of the genetic switch
% Parameters taken from Gardner, Cantor and Collins, Nature, 2000

function dydt = geneticswitch(t, y)

cI = y(1);
lacI = y(2);

% set the max transcription rate in transcripts per second
k_transcription_cI = 0.5;
k_transcription_lacI = 0.5;

% set the leakage transcription rate (ie transcription rate if
% promoter region bound by repressor) in transcripts per second
k_transcription_leakage = 5e-4;

% Set the mRNA and protein degradation rates (per second)
mRNA_half_life = 120; % in seconds
k_mRNA_degradation = log(2)/mRNA_half_life;

% proteins per transcript lifespan
translation_efficiency = 20; 

average_mRNA_lifespan = 1/k_mRNA_degradation; % in seconds

% proteins per transcript per sec
k_translation = translation_efficiency/average_mRNA_lifespan;

% set the Hill coefficients of the repressors
n_cI = 2;
n_lacI = 2;

% Set the dissociation constant for the repressors to their target promoters
% in per molecule per second
KM_tetR = 40;
KM_cI = 40;
KM_lacI = 40;

dydt(1,:) = k_translation * (...
  k_transcription_cI/(1 + (lacI / KM_lacI)^n_lacI) - ...
  k_mRNA_degradation*cI + k_transcription_leakage ...
);

dydt(2,:) = k_translation * (...
  k_transcription_lacI/(1 + (cI / KM_cI)^n_cI) - ...
  k_mRNA_degradation*lacI + k_transcription_leakage ...
);

% Print out some constants for use in lecture
K = k_translation*k_mRNA_degradation
gamma = k_translation * k_transcription_lacI
Kb = (1/KM_lacI)^n_lacI
alpha = gamma * (Kb^(1/n_lacI)) / K