In [1]:
import numpy as np
from scipy.integrate import odeint
import math
import random
import time
import matplotlib
matplotlib.rcParams.update({'font.size': 25})
matplotlib.rc('xtick', labelsize=20)
matplotlib.rc('ytick', labelsize=20)
import matplotlib.pyplot as plt
In [2]:
d_m=math.log(2)/5
d_p=math.log(2)/600
factor=60 #conv sec to min
features={
#RNAP binding free energy, kcal/mol
'dG_RNAP':{"offtarget":-4,
"Phlf":-4.7,
"Srpr":-4.4,
"BetI":-4.6,#
"HlyllR_1":-4.7,
"HlyllR_2":-4.7,
"AmtR":-4.8,
"LacI":-4.5,#
"TetR":-4.5,#
"AraC_1":-4.5,#
"AraC_2":-4.5,#
"Con_A":-4,#
"Con_I":-4,#
"KanR":-4,#
"StrepR":-4#
},#tbd
#TF binding free energy, kcal/mol
'dG_TF':{"offtarget":-6,
"Phlf":-9.2,
"Srpr":-7.4,
"BetI":-8,#
"HlyllR_1":-6.4,
"HlyllR_2":-6.4,
"AmtR":-8.9,
"TetR":-8,#
"AraC_1":-8,#
"AraC_2":-8,#
"LacI":-8#
}, #tbd
#inactive binding free energy
'dG_Inactive':{'TetR':-6.0,#
'LacI':-6.0,#
'AraC':-6.0},#
#mRNA deg rate, 1/min
'd_mRNA':{"Phlf":d_m,#
"Srpr":d_m,#
"BetI":d_m,#
"HlyllR":d_m,#
"AmtR":d_m,#
"YFP":d_m,#
"LacI_TetR":d_m,#
"AraC":d_m,#
"KanR":d_m,#
"StrepR":d_m #
}, #tbd github/cidarlab?
#Translation rate, protein/mRNA/min
'beta':{"Phlf":10,#
"Srpr":10,#
"BetI":10,#
"HlyllR":10,#
"AmtR":10,#
"YFP":10,#
"LacI":10,#
"AraC":10,#
"TetR":10,#
"KanR":10,#
"StrepR":10#
}, #tbd
#number of TF binding sites in promoter
'N_ontarget':{"Phlf":2,#
"Srpr":2,#
"BetI":2,#
"HlyllR_1":2,#
"HlyllR_2":2,#
"AmtR":2,#
"LacI":2,#
"TetR":2,#
"AraC_1":2,#
"AraC_2":2#
}, #tbd
# number of free RNAP
'N_RNAP':750,#
# number of free Ribosomes
'N_RIBO':1000, #
#protein deg rates, 1/min
'd_Protein':{"Phlf":d_p,#
"Srpr":d_p,#
"BetI":d_p,#
"HlyllR":d_p,#
"AmtR":d_p,#
"YFP":d_p,#
"LacI":d_p,#
"AraC":d_p,#
"TetR":d_p,#
"KanR":d_p,#
"StrepR":d_p#
}, #tbd
#growth rate, 1/min
'mu':math.log(2)/45, #tbd
#max promoter transcription rate mRNA/promoter/min
'a_max':{"Phlf":2.84*factor,
"Srpr":2.92*factor,
"BetI":2*factor,#
"HlyllR_1":2.41*factor,
"HlyllR_2":2.41*factor,
"AmtR":1.42*factor,
"LacI":2*factor,#
"TetR":2*factor,#
"AraC_1":2*factor,#
"AraC_2":2*factor,#
"Con_A":2*factor,
"Con_I":2*factor,
"KanR":2*factor,
"StrepR":2*factor
},
#genome length
'genome_len':4600000,
#copy number of circuit
"copy_number":{"Phlf":10,
"Srpr":10,
"BetI":10,
"HlyllR":10,
"AmtR":10,
"YFP":5,
"LacI_TetR":10,
"AraC":10,
"KanR":10,
"StrepR":5},
#cooperativity
'co_op':{'AraC':2,#
'LacI':2,#
'TetR':2#
},
#inducer concentrations used mM
'inducers':{'AraC':0,
'TetR':0,
'LacI':0
},
#kd of inducers
'iKd':{'AraC': 0.005,#
'LacI': 0.05,#
'TetR':0.005 #
},
#kd for hill functions
'Kd':{'AraC':5,#
'LacI':5,#
'TetR':5},#
# temperature
'T':310
}
In [3]:
def Partition_fxn(a_max,dG_RNAP,N_RNAP,dG_TFoff,dG_Active,dG_Inactive,
N_Bindingsites,N_Repressors,fractionActive,genome_len,T):
RT=T*1.987/1000 #temperature(kelvin)*R(kcal/K/mol)
P_RNAP=N_RNAP*math.exp(-dG_RNAP/RT)/genome_len
P_Active=(1 + N_Repressors * fractionActive * math.exp(- (dG_Active-dG_TFoff)) / genome_len)**N_Bindingsites - 1
P_Inactive=(1 + N_Repressors * (1 - fractionActive) * math.exp(-(dG_Inactive-dG_TFoff)) / genome_len)**N_Bindingsites -1;
Pbound=P_RNAP / (1 + P_RNAP + P_Active + P_Inactive)
alpha=a_max*Pbound
return alpha
def Hill_function_R(a_max,Kd,cooperativity,N_Repressor):
alpha = a_max/(1 + (N_Repressor/Kd)**cooperativity);
return alpha
def Hill_function_A(a_max,Ka,cooperativity,N_Activator):
alpha = a_max*((N_Activator)**cooperativity)/(Ka**cooperativity + (N_Repressor)**cooperativity);
return alpha
In [4]:
def rule_30_circuit(y,t,features):
n_strep=y[0]
m_strep=y[1]
p_strep=y[2]
n_kan=y[3]
m_kan=y[4]
p_kan=y[5]
n_ara=y[6]
m_ara=y[7]
p_ara=y[8]
n_lac=y[9]
m_lac=y[10]
p_lac=y[11]
p_tet=y[12]
n_Phlf=y[13]
m_Phlf=y[14]
p_Phlf=y[15]
n_Srpr=y[16]
m_Srpr=y[17]
p_Srpr=y[18]
n_BetI=y[19]
m_BetI=y[20]
p_BetI=y[21]
n_HlyllR=y[22]
m_HlyllR=y[23]
p_HlyllR=y[24]
n_AmtR=y[25]
m_AmtR=y[26]
p_AmtR=y[27]
n_YFP=y[28]
m_YFP=y[29]
p_YFP=y[30]
fractionActive_ara=1/(1+features["inducers"]["AraC"]/features['iKd']['AraC'])
fractionActive_lac=1/(1+features["inducers"]["LacI"]/features['iKd']['LacI'])
fractionActive_tet=1/(1+features["inducers"]["TetR"]/features['iKd']['TetR'])
#print fractionActive_lac
alpha_strep=Partition_fxn(features['a_max']['StrepR'],features['dG_RNAP']['StrepR'],features['N_RNAP'],0,0,0,
0,0,0,features['genome_len'],features['T'])
alpha_kan=Partition_fxn(features['a_max']['KanR'],features['dG_RNAP']['KanR'],features['N_RNAP'],0,0,0,
0,0,0,features['genome_len'],features['T'])
alpha_con_a=Partition_fxn(features['a_max']['Con_A'],features['dG_RNAP']['Con_A'],features['N_RNAP'],0,0,0,
0,0,0,features['genome_len'],features['T'])
alpha_con_L_T=Partition_fxn(features['a_max']['Con_I'],features['dG_RNAP']['Con_I'],features['N_RNAP'],0,0,0,
0,0,0,features['genome_len'],features['T'])
alpha_ara_1=Partition_fxn(features['a_max']['AraC_1'],features['dG_RNAP']['AraC_1'],features['N_RNAP'],features['dG_TF']['offtarget'],features['dG_TF']['AraC_1'],features['dG_Inactive']['AraC'],
features['N_ontarget']['AraC_1'],p_ara,fractionActive_ara,features['genome_len'],features['T'])
alpha_ara_2=Partition_fxn(features['a_max']['AraC_2'],features['dG_RNAP']['AraC_2'],features['N_RNAP'],features['dG_TF']['offtarget'],features['dG_TF']['AraC_2'],features['dG_Inactive']['AraC'],
features['N_ontarget']['AraC_2'],p_ara,fractionActive_ara,features['genome_len'],features['T'])
alpha_tet=Partition_fxn(features['a_max']['TetR'],features['dG_RNAP']['TetR'],features['N_RNAP'],features['dG_TF']['offtarget'],features['dG_TF']['TetR'],features['dG_Inactive']['TetR'],
features['N_ontarget']['TetR'],p_tet,fractionActive_tet,features['genome_len'],features['T'])
alpha_lac=Partition_fxn(features['a_max']['LacI'],features['dG_RNAP']['LacI'],features['N_RNAP'],features['dG_TF']['offtarget'],features['dG_TF']['LacI'],features['dG_Inactive']['LacI'],
features['N_ontarget']['LacI'],p_lac,fractionActive_lac,features['genome_len'],features['T'])
alpha_Phlf=Partition_fxn(features['a_max']['Phlf'],features['dG_RNAP']['Phlf'],features['N_RNAP'],features['dG_TF']['offtarget'],features['dG_TF']['Phlf'],0,
features['N_ontarget']['Phlf'],p_Phlf,1,features['genome_len'],features['T'])
alpha_Srpr=Partition_fxn(features['a_max']['Srpr'],features['dG_RNAP']['Srpr'],features['N_RNAP'],features['dG_TF']['offtarget'],features['dG_TF']['Srpr'],0,
features['N_ontarget']['Srpr'],p_Srpr,1,features['genome_len'],features['T'])
alpha_BetI=Partition_fxn(features['a_max']['BetI'],features['dG_RNAP']['BetI'],features['N_RNAP'],features['dG_TF']['offtarget'],features['dG_TF']['BetI'],0,
features['N_ontarget']['BetI'],p_BetI,1,features['genome_len'],features['T'])
alpha_AmtR=Partition_fxn(features['a_max']['AmtR'],features['dG_RNAP']['AmtR'],features['N_RNAP'],features['dG_TF']['offtarget'],features['dG_TF']['AmtR'],0,
features['N_ontarget']['AmtR'],p_AmtR,1,features['genome_len'],features['T'])
alpha_HlyllR_1=Partition_fxn(features['a_max']['HlyllR_1'],features['dG_RNAP']['HlyllR_1'],features['N_RNAP'],features['dG_TF']['offtarget'],features['dG_TF']['HlyllR_1'],0,
features['N_ontarget']['HlyllR_1'],p_HlyllR,1,features['genome_len'],features['T'])
alpha_HlyllR_2=Partition_fxn(features['a_max']['HlyllR_2'],features['dG_RNAP']['HlyllR_2'],features['N_RNAP'],features['dG_TF']['offtarget'],features['dG_TF']['HlyllR_2'],0,
features['N_ontarget']['HlyllR_2'],p_HlyllR,1,features['genome_len'],features['T'])
#print alpha_lac, alpha_tet ,alpha_ara_1
#print alpha_HlyllR_1
dy=np.zeros(len(y))
dy[0]=0
dy[1]=alpha_strep*n_strep-features['d_mRNA']['StrepR']*m_strep-features['mu']*m_strep
dy[2]=features['beta']['StrepR']*m_strep-features['d_Protein']['StrepR']*p_strep-features['mu']*p_strep
dy[3]=0
dy[4]=alpha_kan*n_kan-features['d_mRNA']['KanR']*m_kan-features['mu']*m_kan
dy[5]=features['beta']['KanR']*m_kan-features['d_Protein']['KanR']*p_kan-features['mu']*p_kan
dy[6]=0
dy[7]=alpha_con_a*n_ara-features['d_mRNA']['AraC']*m_ara-features['mu']*m_ara
dy[8]=features['beta']['AraC']*m_ara-features['d_Protein']['AraC']*p_ara-features['mu']*p_ara
dy[9]=0
dy[10]=alpha_con_L_T*n_lac-features['d_mRNA']['LacI_TetR']*m_lac-features['mu']*m_lac
dy[11]=features['beta']['LacI']*m_lac-features['d_Protein']['LacI']*p_lac-features['mu']*p_lac
dy[12]=features['beta']['TetR']*m_lac-features['d_Protein']['TetR']*p_tet-features['mu']*p_tet
dy[13]=0
dy[14]=(alpha_Srpr+alpha_BetI)*n_Phlf-features['d_mRNA']['Phlf']*m_Phlf-features['mu']*m_Phlf
dy[15]=features['beta']['Phlf']*m_Phlf-features['d_Protein']['Phlf']*p_Phlf-features['mu']*p_Phlf
dy[16]=0
dy[17]=alpha_ara_1*n_Srpr-features['d_mRNA']['Srpr']*m_Srpr-features['mu']*m_Srpr
dy[18]=features['beta']['Srpr']*m_Srpr-features['d_Protein']['Srpr']*p_Srpr-features['mu']*p_Srpr
dy[19]=0
dy[20]=alpha_HlyllR_1*n_BetI-features['d_mRNA']['BetI']*m_BetI-features['mu']*m_BetI
dy[21]=features['beta']['BetI']*m_BetI-features['d_Protein']['BetI']*p_BetI-features['mu']*p_BetI
dy[22]=0
dy[23]=(alpha_tet+alpha_lac)*n_HlyllR-features['d_mRNA']['HlyllR']*m_HlyllR-features['mu']*m_HlyllR
dy[24]=features['beta']['HlyllR']*m_HlyllR-features['d_Protein']['HlyllR']*p_HlyllR-features['mu']*p_HlyllR
dy[25]=0
dy[26]=(alpha_ara_2+alpha_HlyllR_2)*n_AmtR-features['d_mRNA']['AmtR']*m_AmtR-features['mu']*m_AmtR
dy[27]=features['beta']['AmtR']*m_AmtR-features['d_Protein']['AmtR']*p_AmtR-features['mu']*p_AmtR
dy[28]=0
dy[29]=(alpha_Phlf+alpha_AmtR)*n_YFP-features['d_mRNA']['YFP']*m_YFP-features['mu']*m_YFP
dy[30]=features['beta']['YFP']*m_YFP-features['d_Protein']['YFP']*p_YFP-features['mu']*p_YFP
return dy
In [5]:
# Plots the solution to circuit model
def circuit_char(time_points, solution):
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(24, 24))
ax1.set_ylim((0, 3000))
ax2.set_ylim((0, 3000))
ax3.set_ylim((0, 3000))
ax1.set_title('YFP mRNA levels')
ax1.plot(time_points, solution[:, 29], 'orange')
ax1.set_xlabel('Time (minutes)')
ax1.set_ylabel('mRNAs')
ax2.set_title('YFP protein levels')
ax2.plot(time_points, solution[:, 30], 'orange')
ax2.set_xlabel('Time (minutes)')
ax2.set_ylabel('Proteins')
ax3.set_title('TF mRNA levels')
ax3.plot(time_points, solution[:, 14], 'r')
ax3.plot(time_points, solution[:, 17], 'y')
ax3.plot(time_points, solution[:, 20], 'g')
ax3.plot(time_points, solution[:, 23], 'b')
ax3.plot(time_points, solution[:, 26], 'm')
ax3.legend(['Phlf', 'Srpr', 'BetI','HlyllR', 'AmtR'])
ax3.set_xlabel('Time (minutes)')
ax3.set_ylabel('mRNAs')
ax4.set_title('TF protein levels')
ax4.plot(time_points, solution[:, 15], 'r')
ax4.plot(time_points, solution[:, 16], 'y')
ax4.plot(time_points, solution[:, 21], 'g')
ax4.plot(time_points, solution[:, 24], 'b')
ax4.plot(time_points, solution[:, 27], 'm')
ax4.legend(['Phlf', 'Srpr', 'BetI','HlyllR', 'AmtR'])
ax4.set_xlabel('Time (minutes)')
ax4.set_ylabel('Proteins')
fig.subplots_adjust(hspace=0.5)
plt.show()
In [6]:
# Solves circuit model function
def solve_circuit(induction_conditions,features):
t = np.linspace(0, 1000, 1000)
h0 = np.zeros(31)
h0[0]=features['copy_number']['StrepR']
h0[3]=features['copy_number']['KanR']
h0[6]=features['copy_number']['AraC']
h0[9]=features['copy_number']['LacI_TetR']
h0[13]=features['copy_number']['Phlf']
h0[16]=features['copy_number']['Srpr']
h0[19]=features['copy_number']['BetI']
h0[22]=features['copy_number']['HlyllR']
h0[25]=features['copy_number']['AmtR']
h0[28]=features['copy_number']['YFP']
features['inducers']['LacI']=0
features['inducers']['TetR']=0
features['inducers']['AraC']=0
y = odeint(rule_30_circuit, h0, t, (features,))
# circuit_char(t, y)
h0=y[-1,:]
# induction
features['inducers']['LacI']=induction_conditions[0]
features['inducers']['TetR']=induction_conditions[1]
features['inducers']['AraC']=induction_conditions[2]
y = odeint(rule_30_circuit, h0, t, (features,))
circuit_char(t, y)
#solve_circuit([0,0,0],features)
#print "done!"
In [8]:
%matplotlib inline
matplotlib.rcParams.update({'font.size': 25})
solve_circuit([1,0,0],features)
print "done!"
In [ ]:
In [ ]:
In [ ]: