(* Bi 250b - This code is meant to help you calculate the sensitivity to parameters of a constitutive promoter and a negatively-autoregulated promoter *)
(* Unregulated Case *)

(* Define Parameters *)
OpenLoopParams = {
A -> 100 (* mRNA Production Rate *),
G -> 2 (* mRNA Degradation Rate *),
B -> 1000 (* Protein Production Rate *),
D -> 50 (* Protein Degradation Rate *)};

(* expressions for m_dot and p_dot *)
mdot[m_] := (*replace with equation for dm/dt in terms of parameters and m and p *);
pdot[m_,p_] := (*replace with equation for dp/dt in terms of parameters and m and p *);

(* After filling in the expressions for mdot and pdot, hit <Shift+Enter> at the end of the Solve command, this will give you terms for the steady-state values of m and p*) 
Solve[mdot[m]==0 && pdot[m,p] ==0,{m,p}] /. OpenLoopParams

p=(* paste output of p from solve here *);

(* After filling in the expression for p, hit <Shift+Enter> at the end of the NormalizedSensitivityMatrix Line, this will give you the sensitivity for each of the parameters in the order A, G, B and D *) 
SensitivityMatrix = D[p,{{A,G,B,D},1}] /. OpenLoopParams;
NormalizedSensitivityMatrix = 1/p*SensitivityMatrix*{A,G,B,D}/. OpenLoopParams

(* Negative Auto-Regulated Case *)

(* Define Parameters *)
ClosedLoopParams = {
Ar -> 500 (* mRNA Production Rate *),
Gr -> 2 (* mRNA Degradation Rate *),
Br -> 1000 (* Protein Production Rate *),
Dr -> 50 (* Protein Degradation Rate *),
K->  250(* binding-constant of lacI to DNA*)};

(* expressions for m_dot and p_dot *)
mRepressDot[m_,p_] := (*replace with equation for dm/dt in terms of parameters and m and p *);
pRepressDot[m_,p_] := (*replace with equation for dp/dt in terms of parameters and m and p *);

(* After filling in the expressions for mdot and pdot, hit <Shift+Enter> at the end of the Solve command, this will give you terms for the steady-state values of m and p*) 
Solve [mRepressDot[m,p]==0&&pRepressDot[m,p]==0,{m,p}] 
Solve [mRepressDot[m,p]==0&&pRepressDot[m,p]==0,{m,p}] /.ClosedLoopParams

ClosedLoopParams = {
Ar -> 500 (* mRNA Production Rate *),
Gr -> 2 (* mRNA Degradation Rate *),
Br -> 1000 (* Protein Production Rate *),
Dr -> 50 (* Protein Degradation Rate *),
K->  250(* binding-constant of lacI to DNA*),
pRepress ->(* paste numerical output of p from Solve here *)};
pRepress=(* paste symbolic output of p from Solve here *);

(* After filling in the expression for p, hit <Shift+Enter> at the end of the NormalizedSensitivityMatrix Line, this will give you the sensitivity for each of the parameters in the order Ar, Gr, Br, Dr, K and n *) 
SensitivityMatrix = D[pRepress,{{Ar,Gr,Br,Dr, K},1}] /.ClosedLoopParams;
NormalizedSensitivityMatrix = 1/pRepress*SensitivityMatrix*{Ar,Gr,Br,Dr, K}/.ClosedLoopParams