%%
%Start Section 1 - Finding Protein Equilibria Concentrations
V_cytoplasm = 6.7*10^(-16);
V_periplasm = 6.5*10^(-17);

N_A = 6.02215*10^23;

n_p1 = 12;%Copy count for plasmid 1 (copies/cell)
L_p1 = 1400;%Length of pre-mRNA from plasmid 1
n_p2 = 5;%Copy count for plasmid 2 (copies/cell)
p2_conc = n_p2/(V_cytoplasm*N_A);
L_p2 = 3430;

k_transcription = 70; %transcription rate (70 nucleotides/sec)
k_translation  = 40; %translation rate (40 amino acids/sec)
k_mRNA_degradation = 0.0058; %/s (Newcastle)
k_protein_degradation = 0.0012; %/s for MerR, MerT, MerP, MT (Newcastle)

L_MerR = 145; %Length of MerR protein (# of amino acids)
L_MerP = 276; %Length of MerP protein (# of amino acids)
L_MerT = 116; %Length of MerT protein (# of amino acids)
L_MetT = 73; %Length of MetT protein (# of amino acids)
L_Ag43 = 1040; %Length of Ag43 protein (# of amino acids)


[T1 MERR_MRNA]=ode45(@(t, MerR_mRNA)(0.1*k_transcription*n_p1/L_p1-k_mRNA_degradation*MerR_mRNA),0:10:10000,0);
d1=ode45(@(t, MerR_mRNA)(0.1*k_transcription*n_p1/L_p1-k_mRNA_degradation*MerR_mRNA),0:10.000000:10000,0);
plot(T1, MERR_MRNA)
legend ('n(MerR_mRNA)/cell')
hold on;

figure
[T2 MERR]=ode45(@(t, MerR)((k_translation*deval(d1,t)/L_MerR)/(N_A*V_cytoplasm)-k_protein_degradation*MerR), 0:10.000000:10000,0);
d2=ode45(@(t, MerR)((k_translation*deval(d1,t)/L_MerR)/(N_A*V_cytoplasm)-k_protein_degradation*MerR), 0:10.000000:10000,0);
semilogy(T2, MERR,'k')
hold on;
[T3 MERP]=ode45(@(t, MerP)((k_translation*deval(d1,t)/L_MerP)/(N_A*V_periplasm)-k_protein_degradation*MerP), 0:10.000000:10000,0);%note expressed in periplasm
semilogy(T3, MERP,'b')
[T4 MERT]=ode45(@(t, MerT)((k_translation*deval(d1,t)/L_MerT)/(N_A*V_periplasm)-k_protein_degradation*MerT), 0:10.000000:10000,0);%note expressed in periplasm
semilogy(T4, MERT,'m')
[T5 METT]=ode45(@(t, MetT)((k_translation*deval(d1,t)/L_MetT)/(N_A*V_cytoplasm)-k_protein_degradation*MetT), 0:10.000000:10000,0);
semilogy(T5, METT,'g')


kMerR2_p2_A = 500;
kMerR2_p2_D = 0.0005;

[T6 MERR2_P2]=ode45(@(t, MerR2_p2)(kMerR2_p2_A*(0.5*deval(d2,t)-MerR2_p2)*(p2_conc-MerR2_p2)-kMerR2_p2_D*MerR2_p2), 0:10.000000:10000,0);

semilogy(T6,MERR2_P2, 'c')
legend ('[MerR]','[MerP]','[MerT]','[MetT]', '[MerR_2.p2]');


Total_MerR=MERR(1000);
Total_MerP=MERP(1000);
Total_MerT=MERT(1000);
Total_MetT=METT(1000);
Total_MerR2_p2=MERR2_P2(1000);

%End Section 1

%%
%define variables

Hg_total = 100*10^(-6);% Initial Hg concentration (M)
OD = 100; %Optical Density = A600
Cell_density = OD*10^12; %Cell density in cells/L
V_cytoplasm = 6.7*10^(-16);%Volume of cytoplasm of 1 cell (L)
V_periplasm = 6.5*10^(-17);%Volume of periplasm of 1 cell (L)

Kd_HgMerP = 2.7*10^(-6);%Dissociation constant for binding of Hg to MerP (M)
k_fwd_HgMerP = 10;%Guess at speed of reaction (forward)
k_bwd_HgMerP = k_fwd_HgMerP*Kd_HgMerP;%Reverse reaction constant to maintian Kd.

k_uptake = 10000;%Uptake rate constant of Hg through MerP-MerT pathway
k_efflux = 0.01;%Efflux(diffusion) rate constant of Hg through inner membrane

k_rate_HgMerR2p2 = 0.05;%Guess at speed of reaction (based on distance from equilibrium {R = k(X(eq)-X)}
n_H = 2.6; %Hill coefficient for binding of Hg to MerR2_p2 complex.
K_M = 1.2*10^(-8);%Michaelis constant for binding of Hg to MerR2_p2 complex (M)

n_sites = 8;%Estimated number of binding sites on MetT
k_fwd_HgMetT = 10000;%Guess forward rate constant Hg-MetT binding /(M.s)]
k_bwd_HgMetT = 0.001;%Guess reverse Reaction rate Hg-MetT binding [/s]

%%
%ODEs
%Hg_int = y =v(1)
%Hg_ext = x =(Hg_total – Cell_density*(V_cytoplasm*y + V_periplasm*a))
%Hg_free = z = (y-b-n_sites*c)
%Hg_MerR2p2 = b = v(2)
%Hgn_MetT = c = v(3)
%Hg_MerP = a = v(4)
%Ag43_mRNA = e = v(5)
%Ag43 = q = v(6)

%Hg_int= dy/dt = f(x,y,z,a,b,c,e,q)
%Hg_MerR2p2 = db/dt = g(x,y,z,a,b,c,e,q)
%Hgn_MetT = dc/dt = h(x,y,z,a,b,c,e,q)
%Hg_MerP = da/dt = j(x,y,z,a,b,c,e,q)
%Ag43_mRNA = de/dt = l(x,y,z,a,b,c,e,q)
%Ag43 = dq/dt = m(x,y,z,a,b,c,e,q)

%f = (k_uptake*v(4)*Total_MerT-k_efflux*(v(1)-v(2)-n_sites*v(3)));
%g = (k_rate_HgMerR2p2*((Total_MerR2_p2*((v(1)-v(2)-n_sites*v(3))^n_H)/( (v(1)-v(2)-n_sites*v(3))^(n_H)+K_M^(n_H)))-v(2)));
%h = (k_fwd_HgMetT*(v(1)-v(2)-n_sites*v(3))*(Total_MetT-v(3))-k_bwd_HgMetT*v(3));
%j = (k_fwd_HgMerP*(Hg_total-Cell_density*(V_cytoplasm*v(1) + V_periplasm*v(4)))*(Total_MerP-v(4))-k_bwd_HgMerP*v(4));
%l = (k_transcription*v(2)/L_p2-k_mRNA_degradation*v(5));
%q = (k_translation*v(5)/L_Ag43-k_protein_degradation*v(6));

[T V]=ode45(@(t,v)([(k_uptake*v(4)*Total_MerT-k_efflux*(v(1)-v(2)-n_sites*v(3))); (k_rate_HgMerR2p2*((Total_MerR2_p2*((v(1)-v(2)-n_sites*v(3))^n_H)/( (v(1)-v(2)-n_sites*v(3))^(n_H)+K_M^(n_H)))-v(2))); (k_fwd_HgMetT*(v(1)-v(2)-n_sites*v(3))*(Total_MetT-v(3))-k_bwd_HgMetT*v(3)); (k_fwd_HgMerP*(Hg_total-Cell_density*(V_cytoplasm*v(1) + V_periplasm*v(4)))*(Total_MerP-v(4))-k_bwd_HgMerP*v(4)); (k_transcription*v(2)/L_p2-k_mRNA_degradation*v(5)); (k_translation*v(5)/L_Ag43-k_protein_degradation*v(6));]), 0:1:10000, [0, 0, 0, 0, 0, 0]);

T7=0:1:10000;
i=1;
while i<=10001
    Hg_ext(i)=(Hg_total-Cell_density*(V_cytoplasm*V(i,1)+V_periplasm*V(i,4)));
    i=i+1;
end

i=1;
while i<=10001
    Hg_free(i)=(V(i,1)-V(i,2)-n_sites*V(3));
    i=i+1;
end

figure
semilogy(T,V(:,1),'k')%Hg_int
hold on;
semilogy(T,V(:,2),'b')%Hg_MerR2_p2
semilogy(T,V(:,3),'g')%Hgn_MetT
semilogy(T,V(:,4),'m')%Hg_MerP
semilogy(T7,Hg_ext,'y')%Hg_ext
semilogy(T7,Hg_free,'r')%Hg_free
semilogy(T,V(:,5),'c')%Ag43_mRNA
semilogy(T,V(:,6),'w')%Ag43
legend('[Hg]_i_n_t)','[Hg.MerR_2.p2]','[Hg_n.MetT]','[Hg.MerP]','[Hg]_e_x_t','[Hg]_f_r_e_e','[Ag43.mRNA]','[Ag43]');