In [1]:
%pylab inline
import tellurium as te
from SloppyCell.ReactionNetworks import Utility
from scipy.stats import linregress

Populating the interactive namespace from numpy and matplotlib
In [2]:
#For reading excel files used in the modelling
!pip install xlwt

Requirement already satisfied: xlwt in /usr/local/lib/python2.7/dist-packages
In [3]:
import os
from pandas import *
import matplotlib.pyplot
from scipy.integrate import odeint
from lmfit import minimize, Parameter, Parameters, report_fit
from scipy.interpolate import interp1d, CubicSpline

##This allows is for writing in excel the results

from tempfile import TemporaryFile
from xlwt import Workbook

book = Workbook()

##########################

%run simple_model.py
%run parameters.py
In [4]:
def residuals(params, data, inter):
    y = array([params['y']])
    prediction=odeint(simple_model,y,range(0,24*7+2,1),args=(u1,params['degradation'], inter))
    p = {}
    for idx, i in enumerate(range(0,24*7+2,1)):
        p[i] = prediction[idx][0]

    res = []
    B_num = 0
    B_den = 0
    for k  in data.keys():
        B_num+=(p[k]*data[k])
        B_den+= p[k]**2
    B = B_num/B_den



    for k  in data.keys():
        res.append((B*p[k]-data[k])**2)
    
    return res

def sf(params, data, inter):
    y = array([params['y']])
    prediction=odeint(simple_model,y,range(0,24*7+2,1),args=(u1,params['degradation'], inter))
    p = {}
    for idx, i in enumerate(range(0,24*7+2,1)):
        p[i] = prediction[idx][0]

    res = []
    B_num = 0
    B_den = 0
    for k  in data.keys():
        B_num+=(p[k]*data[k])
        B_den+= p[k]**2
    B = B_num/B_den



    for k  in data.keys():
        res.append((B*p[k]-data[k])**2)
    
    return B

The RNA data is loaded first. I use the standard that I developed for my PhD. First we have the RNA data loaded into sloppy cell

In [5]:
book = ExcelFile('PhotoPeriodTiMetEdited.xlsx')
book.sheet_names

data = {}
for network in book.sheet_names:
    data[network] =book.parse(network)

genes = {'cL_m':'lhy','cC_m':'cca1','cP9_m':'prr9','cP7_m':'prr7','cP5_m':'prr5',
         'cT_m':'toc1','cLUX_m':'lux','cG_m':'gi','cE3_m':'elf3','cE4_m':'elf4'}

genotypes = {}

##From here we only use a subset of the data we are interested in constant light

#genotypes['pp6-18'] = []
#genotypes['pp8-16'] = []
#genotypes['pp12-12'] = []
#genotypes['pp18-6'] = []
genotypes['col_0_LL'] = []
#genotypes['col_0_DD'] = []
genotypes['prr79_LL'] = ['cP9_m','cP7_m']
#genotypes['prr79_DD'] = ['cP9_m','cP7_m']
genotypes['lhycca1_LL'] = ['cL_m']
#genotypes['lhycca1_DD'] = ['cL_m']
genotypes['toc1_LL'] = ['cT_m']
#genotypes['toc1_DD'] = ['cT_m']
genotypes['gi_LL'] = ['cG_m']
#genotypes['gi_DD'] = ['cG_m']


nuclei_number_gFW = 25.0e6
uncertainty = 0.6
displacement=24*10


## This function will be introduced in SloppyCell Utility module cos it is a very in order to make script more compact

def data_formater(data,genotypes,genes ,nuclei_number_gFW, displacement):
    '''
    Formating of data for sloppycell starting from an excel file
    This formating assumes that it is being fed with a dictionary of panda Data frames

    '''

    network_data={}

    for network in data.keys():
        if network not in genotypes.keys():
            continue

        gene_dict = {}
        for gene in genes.keys():
            try:
                if gene in genotypes[network]:
                    continue
            except:
                pass
            time_dict={}
            for zt in range(0,len(data[network][genes[gene]]),2):

                if not (np.isnan(data[network][genes[gene]][zt]) or np.isnan(data[network][genes[gene]][zt+1])):
                    points  = np.log(np.array([data[network][genes[gene]][zt],
                                               data[network][genes[gene]][zt+1]])/nuclei_number_gFW)
                    data_mean=np.mean(points)
                    data_var=np.std(points)                    
                else:
                    if np.isnan(data[network][genes[gene]][zt+1]) and not np.isnan(data[network][genes[gene]][zt]):
                        points  = np.log(np.array([data[network][genes[gene]][zt]])/nuclei_number_gFW)
                        data_mean=np.mean(points)
                        data_var = 0.6
                    elif np.isnan(data[network][genes[gene]][zt]) and not np.isnan(data[network][genes[gene]][zt+1]):
                        points  = np.log(np.array([data[network][genes[gene]][zt+1]])/nuclei_number_gFW)
                        data_mean=np.mean(points)
                        data_var = 0.6



                time_dict[zt+displacement] = (data_mean,data_var)
            gene_dict['log_'+gene]=time_dict
        network_data[network]=gene_dict
    return network_data

network_data = data_formater(data, genotypes,genes,nuclei_number_gFW, displacement)
In [6]:
exp(network_data['col_0_LL']['log_cL_m'][240][0])/exp(network_data['col_0_LL']['log_cC_m'][240][0])
Out[6]:
3.3138024885658273
In [7]:
figure(dpi=600)
def interpol_data(network_data, network, var):
    temp_time = array(sorted(network_data[network][var].keys()))
    temp = []
    for i in temp_time:
        temp.append(network_data[network][var][i][0])
    temp = array(temp)
    return temp_time-24*10-24, temp

interpolations = {}
for i in network_data['col_0_LL'].keys():
    temp_time, temp = interpol_data(network_data, 'col_0_LL', i)
    interpolations[i] = CubicSpline(temp_time,temp)

selected_genes = ['cL_m', 'cC_m']
colours = ['red','blue','orange','green']
k = 0
for j in selected_genes:
    plot(linspace(0,24*3,2000)-24,(interpolations['log_'+j](linspace(0,24*3,2000)-24)), lw=2,color=colours[k], label=genes[j].upper()+' mRNA')
    for i in array(sorted(network_data['col_0_LL']['log_'+j].keys())):
        errorbar(i-24*10-24, network_data['col_0_LL']['log_'+j][i][0], yerr=network_data['col_0_LL']['log_'+j][i][1], lw=2,fmt='s', color=colours[k])
    k+=1
ylim(-4,9)
xlim(-24,48)
xticks(range(0,48,12), fontsize=16)
yticks(fontsize=16)
ylabel('ln(#transcripts/cell)', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
legend(loc='lower right', fontsize=16)
axvspan(-12,0, color='gray', alpha=0.7)
for i in range(0,3,1):
    axvspan(12+24*i,24*(i+1), color='gray', alpha=0.1)

tight_layout()
#savefig('cubic_interpolation_example.svg', format='svg', dpi=600, transparent=True)

Interpolation of protein CCA1 in entrainement conditions

In [8]:
figure(dpi=600)
temp_time, temp = interpol_data(network_data, 'col_0_LL', 'log_cL_m')
b = tile(temp[0:12],5)
b = concatenate((b,temp[12:]))
lhyinter = CubicSpline(range(0,len(b)*2,2),b)
plot(range(0,len(b)*2,2),exp(b))
plot(linspace(2,24*7,200),exp(lhyinter(linspace(2,24*7,200))))
xlim(24*3,24*7)
xlabel('Time in constant light (h)', fontsize=22)
axvspan(-12,0, color='gray', alpha=0.7)
for i in range(0,3,1):
    axvspan(12+24*i,24*(i+1), color='gray', alpha=0.1)
In [9]:
temp_time, temp = interpol_data(network_data, 'col_0_LL', 'log_cC_m')
b = tile(temp[0:12],5)
b = concatenate((b,temp[12:]))
cca1inter = CubicSpline(range(0,len(b)*2,2),b)
plot(range(0,len(b)*2,2),exp(b))
plot(linspace(2,24*7,200),exp(cca1inter(linspace(2,24*7,200))))
Out[9]:
[<matplotlib.lines.Line2D at 0x407b2161d0>]

Simple model that has one differential equation. We do not fitt the model to the data as there are degradation rates for LHY. The degradation rate comes from Hansen Louise un published

In [10]:
from scipy.integrate import odeint
def simple_model(y,t,u1,d1,mdata1):
    mRNA = mdata1(t)
    Prot0 = y[0]
    dydt=range(1)
    dydt[0]= u1*exp(mRNA)-d1*Prot0
    return dydt
In [11]:
#in aminoacids 
amin_sec = 3.
rib_tranc = 10.
CCA1_length = 609.
CCA1_trans_rate = 1/(CCA1_length/amin_sec/10)*3600
print CCA1_trans_rate
amin_sec = 3.
rib_tranc = 10.
LHY_length = 646.
LHY_trans_rate = 1/(LHY_length/amin_sec/10)*3600
print LHY_trans_rate

177.339901478
167.182662539

Model behaviour in 12L:12D using the last LD cycle of Wang figure 1

In [12]:
figure(dpi=600)
deg_cca1 = read_excel('CCA1_degradation_louise.xlsx', sheet_name='Sheet1')
CCA1_reg = linregress(deg_cca1.iloc[:,0]/60.,log(mean(deg_cca1.iloc[:,1:4], axis=1)/mean(deg_cca1.iloc[:,1:4], axis=1)[0]))
plot(deg_cca1.iloc[:,0]/60.,log(mean(deg_cca1.iloc[:,1:4], axis=1)/mean(deg_cca1.iloc[:,1:4], axis=1)[0]))
plot(deg_cca1.iloc[:,0]/60., CCA1_reg.slope*deg_cca1.iloc[:,0]/60.+CCA1_reg.intercept)
xlabel('Time in hours', fontsize=16)
ylabel('ln(CCA1 signal/(CCA1 t0))', fontsize=16)
Out[12]:
Text(0,0.5,'ln(CCA1 signal/(CCA1 t0))')
In [14]:
deg_cca1 = read_excel('CCA1_degradation_louise.xlsx', sheet_name='Sheet2')
lhy_reg = linregress(deg_cca1.time/60,deg_cca1.deg)
lhy_reg.slope
Out[14]:
-1.7688
In [15]:
pCCA1 = read_excel('ProteinTimeseriesDB2.xlsx', sheet_name='pCCA1(Wang1998)')
In [16]:
res = {}
pro_data = {}
In [17]:
print log(2)/lhy_reg.slope*(-1)
print log(2)/CCA1_reg.slope*(-1)

0.3918742540479112
0.20434558391843624
In [18]:
figure(dpi=600)
y = zeros(1)
t = range(0,168,1)
t2 = numpy.array(range(1,168,1))
res['LHY'] = odeint(simple_model,y,t,args=(LHY_trans_rate,lhy_reg.slope*(-1),lhyinter))
#plot(pCCA1.time,log(pCCA1.pCCA1*10000), 'bo-', label='Knowles2008')

xticks(range(0,24*4,12))
plot(t2-24*5,res['LHY'][1:,0], lw=3, color= 'red', label='pLHY')
yticks(fontsize=16)
xticks(fontsize=16)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in h', fontsize=22)
xticks(fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(24*2+12,24*3, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
y = zeros(1)
res['CCA1'] = odeint(simple_model,y,t,args=(CCA1_trans_rate,CCA1_reg.slope*(-1),cca1inter))
plot(t2-24*5,(res['CCA1'][1:,0]), lw=3, color= 'blue', label='pCCA1')
xlim(-24,24*2-4)
yticks(fontsize=16)
xticks(fontsize=16)
ylim(1,1e6)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
xticks(fontsize=16)
plot(t2-24*5,(res['LHY'][1:,0]+res['CCA1'][1:,0]), lw=3, color= 'orange', label='pLHY+pCCA1')

pro_data['cL'] = {}
for idx, i in enumerate(res['CCA1']):
    if idx >= 24*4:
        pro_data['cL'][t[idx]] =  (res['CCA1'][idx][0]+res['LHY'][idx][0],0.1)

yscale('log')
legend(loc='upper right', fontsize=16)
tight_layout()
#savefig('lhy_cca1_simple_model.svg', format='svg', dpi=600, transparent=True)
In [19]:
def simple_model_ligh_deg(y,t,u1,d1,d2,period,dusk,tw,mdata1):
    if t < 24*5:
        l=0.5*((1+np.tanh((t-period*np.floor(t/period))/tw))-(1+np.tanh((t-period*np.floor(t/period)-dusk)/tw))+(1+np.tanh((t-period*np.floor(t/period)-period)/tw)))
    else :
        l=1.0
    mRNA = mdata1(t)
    Prot0 = y[0]
    dydt=range(1)
    dydt[0]= u1*exp(mRNA)-(d1+d2*(1-l))*Prot0
    return dydt
In [20]:
prr7decay = read_excel('data/prr7prr9_decay.xlsx',sheet_name='PRR7decaycurve')
     
PRR7_reg_light = linregress(prr7decay.timeH,log(prr7decay.ZT12L))
plot(prr7decay.timeH,log(prr7decay.ZT12L))
plot(prr7decay.timeH, PRR7_reg_light.slope*prr7decay.timeH+PRR7_reg_light.intercept)

PRR7_reg_dark = linregress(prr7decay.timeH,log(prr7decay.ZT12D))
plot(prr7decay.timeH,log(prr7decay.ZT12D))
plot(prr7decay.timeH, PRR7_reg_dark.slope*prr7decay.timeH+PRR7_reg_dark.intercept)
Out[20]:
[<matplotlib.lines.Line2D at 0x40995b0390>]
In [21]:
temp_time, temp = interpol_data(network_data, 'col_0_LL', 'log_cP7_m')
b = tile(temp[0:12],5)
b = concatenate((b,temp[12:]))
p7inter = CubicSpline(range(0,len(b)*2,2),b)
plot(range(0,len(b)*2,2),exp(b))
plot(linspace(2,24*7,200),exp(p7inter(linspace(2,24*7,200))))
Out[21]:
[<matplotlib.lines.Line2D at 0x407ce20210>]
In [22]:
figure(dpi=600)
amin_sec = 3.
rib_tranc = 10.
PRR7_length = 727.
PRR7_trans_rate = 1/(PRR7_length/amin_sec/10)*3600
print PRR7_trans_rate
y = zeros(1)
period = 24.
dusk = 12.
tw = 1.
daypat = 'LD'
t = numpy.array(range(0,24*7,1))
res['PRR7'] = odeint(simple_model_ligh_deg,y,t,
              args=(PRR7_trans_rate,PRR7_reg_light.slope*(-1),PRR7_reg_dark.slope*(-1)-PRR7_reg_light.slope*(-1),
                    period,dusk,tw,
                    p7inter))
pro_data['cP7'] = {}
for idx, i in enumerate(res['PRR7']):
    if idx >= 24*4:
        pro_data['cP7'][t[idx]] =  (res['PRR7'][idx][0],0.1)

        
plot(t-24*5,res['PRR7'][:,0], lw=3, color= 'black', label='LHY')
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in h', fontsize=22)
xticks(fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(12,24, color='gray',alpha=0.1)
axvspan(24+12,24*3, color='gray',alpha=0.1)
xticks(range(0,24*3,12),fontsize=16)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in contant light (h)', fontsize=22)
xlim(-24,24*2)
yticks(fontsize=16)

148.555708391
Out[22]:
(array([-10000.,      0.,  10000.,  20000.,  30000.,  40000.,  50000.,
         60000.,  70000.]), <a list of 9 Text yticklabel objects>)
In [23]:
print (pro_data['cP7'][96+9][0]-pro_data['cP7'][96+24][0])/2
print (pro_data['cP7'][96+9+24][0]-pro_data['cP7'][96+24*2+1][0])/2
print (pro_data['cP7'][96+9+24][0]-pro_data['cP7'][96+24*2+1][0])/(pro_data['cP7'][96+9][0]-pro_data['cP7'][96+24][0])

15411.987765454134
24538.23309287341
1.5921523859417859
In [24]:
prr7decay = read_excel('data/prr7prr9_decay.xlsx',sheet_name='PRR9decay')
     
PRR7_reg_light = linregress(prr7decay.timemin,log(prr7decay.morning))
plot(prr7decay.timemin,log(prr7decay.morning))
plot(prr7decay.timemin, PRR7_reg_light.slope*prr7decay.timemin+PRR7_reg_light.intercept)
Out[24]:
[<matplotlib.lines.Line2D at 0x409c4a8c50>]
In [26]:
temp_time, temp = interpol_data(network_data, 'col_0_LL', 'log_cP9_m')
b = tile(temp[0:12],5)
b = concatenate((b,temp[12:]))
p7inter = CubicSpline(range(0,len(b)*2,2),b)
plot(range(0,len(b)*2,2),exp(b))
plot(linspace(2,24*7,200),exp(p7inter(linspace(2,24*7,200))))
Out[26]:
[<matplotlib.lines.Line2D at 0x40a06a4650>]
In [27]:
figure(figsize=(10,7))
amin_sec = 3.
rib_tranc = 10.
PRR7_length = 468.
PRR7_trans_rate = 1/(PRR7_length/amin_sec/10)*3600
print PRR7_trans_rate
y = zeros(1)
period = 24.
dusk = 12.
tw = 1.
daypat = 'LD'
t = numpy.array(range(0,24*7,1))
res['PRR9'] = odeint(simple_model,y,t,
              args=(PRR7_trans_rate,PRR7_reg_light.slope*(-1),
                    p7inter))
plot(t-24*5,res['PRR9'][:,0], lw=3, color= 'black', label='LHY')
pro_data['cP9'] = {}
for idx, i in enumerate(res['PRR9']):
    if idx >= 24*4:
        pro_data['cP9'][t[idx]] =  (res['PRR9'][idx][0],0.1)

xticks(fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.1)
axvspan(24*2+12,24*3, color='gray',alpha=0.1)
axvspan(12,24, color='gray',alpha=0.1)
xticks(range(0,24*3,12),fontsize=16)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in contant light (h)', fontsize=22)
xlim(-24,24*2)
yticks(fontsize=16)
axvspan(12,24, color='gray',alpha=0.7)
axvline(2)
axvline(8)

230.769230769
Out[27]:
<matplotlib.lines.Line2D at 0x40a0a208d0>
In [28]:
temp_time, temp = interpol_data(network_data, 'col_0_LL', 'log_cP5_m')
b = tile(temp[0:12],5)
b = concatenate((b,temp[12:]))
p7inter = CubicSpline(range(0,len(b)*2,2),b)
plot(range(0,len(b)*2,2),exp(b))
plot(linspace(2,24*7,200),exp(p7inter(linspace(2,24*7,200))))

prr7decay = read_excel('data/prr7prr9_decay.xlsx',sheetname='PRR5decay')
     
PRR7_reg_wl = linregress(prr7decay.time,log(prr7decay.wl))
plot(prr7decay.time,log(prr7decay.wl))
plot(prr7decay.time, PRR7_reg_wl.slope*prr7decay.time+PRR7_reg_wl.intercept)

PRR7_reg_dark = linregress(prr7decay.time,log(prr7decay.dark))
plot(prr7decay.time,log(prr7decay.dark))
plot(prr7decay.time, PRR7_reg_dark.slope*prr7decay.time+PRR7_reg_dark.intercept)

figure(figsize=(10,7))
amin_sec = 3.
rib_tranc = 10.
PRR7_length = 558.
PRR7_trans_rate = 1/(PRR7_length/amin_sec/10)*3600
print PRR7_trans_rate
y = zeros(1)
period = 24.
dusk = 12.
tw = 1.
daypat = 'LD'
t = numpy.array(range(0,24*7,1))
res['PRR5'] = odeint(simple_model_ligh_deg,y,t,
              args=(PRR7_trans_rate,PRR7_reg_wl.slope*(-1),PRR7_reg_dark.slope*(-1)-PRR7_reg_wl.slope*(-1),
                    period,
                    dusk,
                    tw,
                    p7inter))

plot(t-24*5,res['PRR5'][:,0], lw=3, color= 'black', label='LHY')
pro_data['cP5'] = {}
for idx, i in enumerate(res['PRR5']):
    if idx >= 24*4:
        pro_data['cP5'][t[idx]] =  (res['PRR5'][idx][0],0.1)

xticks(fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.1)
axvspan(24*2+12,24*3, color='gray',alpha=0.1)
axvspan(12,24, color='gray',alpha=0.1)
xticks(range(0,24*3,12),fontsize=16)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in contant light (h)', fontsize=22)
xlim(-24,24*2)
yticks(fontsize=16)

/usr/local/lib/python2.7/dist-packages/pandas/util/_decorators.py:188: FutureWarning:

The `sheetname` keyword is deprecated, use `sheet_name` instead


193.548387097
Out[28]:
(array([-10000.,      0.,  10000.,  20000.,  30000.,  40000.,  50000.,
         60000.,  70000.]), <a list of 9 Text yticklabel objects>)
In [29]:
temp_time, temp = interpol_data(network_data, 'col_0_LL', 'log_cT_m')
b = tile(temp[0:12],5)
b = concatenate((b,temp[12:]))
p7inter = CubicSpline(range(0,len(b)*2,2),b)
plot(range(0,len(b)*2,2),exp(b))
plot(linspace(2,24*7,200),exp(p7inter(linspace(2,24*7,200))))
figure(figsize=(10,7))
prr7decay = read_excel('data/prr7prr9_decay.xlsx',sheetname='TOC1decay')
     
PRR7_reg_wl = linregress(prr7decay.time,log(prr7decay.normdecay))
plot(prr7decay.time,log(prr7decay.normdecay))
plot(prr7decay.time, PRR7_reg_wl.slope*prr7decay.time+PRR7_reg_wl.intercept)

figure(figsize=(10,7))
amin_sec = 3.
rib_tranc = 10.
PRR7_length = 618.
PRR7_trans_rate = 1/(PRR7_length/amin_sec/10)*3600
print PRR7_trans_rate
y = zeros(1)
period = 24.
dusk = 12.
tw = 1.
daypat = 'LD'
t = numpy.array(range(0,24*7,1))
res['TOC1'] = odeint(simple_model,y,t,
              args=(PRR7_trans_rate,PRR7_reg_wl.slope*(-1),
                    p7inter))
plot(t-24*5,res['TOC1'][:,0], lw=3, color= 'black', label='TOC1')

pro_data['cT'] = {}
for idx, i in enumerate(res['TOC1']):
    if idx >= 24*4:
        pro_data['cT'][t[idx]] =  (res['TOC1'][idx][0],0.1)

xticks(fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.1)
axvspan(24*2+12,24*3, color='gray',alpha=0.1)
axvspan(12,24, color='gray',alpha=0.1)
xticks(range(0,24*3,12),fontsize=16)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in contant light (h)', fontsize=22)
xlim(-24,24*2)
yticks(fontsize=16)

174.757281553
Out[29]:
(array([-5000.,     0.,  5000., 10000., 15000., 20000., 25000., 30000.,
        35000.]), <a list of 9 Text yticklabel objects>)
In [31]:
!ls data

CCA1WangLD1998		       PBM_vs_U20172.xlsx
CCA1_degradation_louise.xlsx   PhotoPeriodTiMetEdited.xlsx
ELF3ELF4LUXProtein.xlsx        UrquizaGarcia_ProteinTimeseriesDB2.xlsx
LUX_affinities_PBM_U2017.xlsx  prr7prr9_decay.xlsx
In [30]:
figure(dpi=600)
t = numpy.array(range(0,24*7,1))
plot(t-24*5,res['PRR9'][:,0], lw=3, color= 'black', label='P9')
pProtein = read_excel('data/UrquizaGarcia_ProteinTimeseriesDB2.xlsx', sheetname='PRR9(Fujiware2008fig1F)')
plot(pProtein.time,pProtein.pPRR9*9000, 'ko-')
w_data={}
w_data['cP9'] ={}
for idx,i in enumerate(pProtein.time):
     w_data['cP9'][int(pProtein.iloc[idx,0])+24*5]  =  pProtein.iloc[idx,1]

plot(t-24*5,res['PRR7'][:,0], lw=3, color= 'blue', label='P7')
pProtein = read_excel('data/UrquizaGarcia_ProteinTimeseriesDB2.xlsx', sheetname='pPRR7(fujiware2008fig1F)')
plot(pProtein.time,pProtein.pPRR7*60000, 'bo-')
w_data['cP7'] ={}
for idx,i in enumerate(pProtein.time):
     w_data['cP7'][int(pProtein.iloc[idx,0])+24*5]  =  pProtein.iloc[idx,1]
xticks(fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(24*2+12,24*3, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
xticks(range(0,24*3,12),fontsize=16)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in contant light (h)', fontsize=22)
xlim(-24,24*2)
yticks(fontsize=16)
ylim(1,1e6)
yscale('log')
legend(loc='lower left', fontsize=16)
tight_layout()
#savefig('PRR9_PRR7_protein_prediction.svg',format='svg',dpi=600,transparent=True)
#yscale('log')
In [31]:
figure(dpi=600)
plot(t-24*5,res['PRR5'][:,0], lw=3, color= 'orange', label='P5')
pProtein = read_excel('data/UrquizaGarcia_ProteinTimeseriesDB2.xlsx', sheetname='pPRR5(Fujiware2008fig1F)')
plot(pProtein.time,pProtein.pPRR5*60000, 'o-', color='orange')
w_data['cP5'] ={}
for idx,i in enumerate(pProtein.time):
     w_data['cP5'][int(pProtein.iloc[idx,0])+24*5]  =  pProtein.iloc[idx,1]

plot(t-24*5,res['TOC1'][:,0], lw=3, color= 'red', label='TOC1')
pProtein = read_excel('data/UrquizaGarcia_ProteinTimeseriesDB2.xlsx', sheetname='pTOC1(Fujiware2008fig1F)')
plot(pProtein.time,pProtein.pTOC1*30000, 'ro-')
w_data['cT'] ={}
for idx,i in enumerate(pProtein.time):
     w_data['cT'][int(pProtein.iloc[idx,0])+24*5]  =  pProtein.iloc[idx,1]
xticks(fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(24*2+12,24*3, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
xticks(range(0,24*3,12),fontsize=16)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in contant light (h)', fontsize=22)
xlim(-24,24*2)
ylim(1,1e6)
yticks(fontsize=16)
legend(loc='lower right', fontsize=16)
yscale('log')
tight_layout()
#savefig('PRR5_TOC1_protein_prediction.svg',format='svg',dpi=600,transparent=True)

#savefig('figures/PRR51_simple_model.png', format='png', dpi=300)
In [32]:
def simple_model(y,t,u1,d1,mdata1):
    mRNA = mdata1(t)
    Prot0 = y[0]
    dydt=range(1)
    dydt[0]= u1*exp(mRNA)-d1*Prot0
    return dydt
In [34]:
def residuals(params, data, inter):
    y = array([params['y']])
    prediction=odeint(simple_model,y,range(0,24*7+2,1),args=(u1,params['degradation'], inter))
    p = {}
    for idx, i in enumerate(range(0,24*7+2,1)):
        p[i] = prediction[idx][0]

    res = []
    B_num = 0
    B_den = 0
    for k  in data.keys():
        B_num+=(p[k]*data[k])
        B_den+= p[k]**2
    B = B_num/B_den



    for k  in data.keys():
        res.append((B*p[k]-data[k])**2)
    
    return res

def sf(params, data, inter):
    y = array([params['y']])
    prediction=odeint(simple_model,y,range(0,24*7+2,1),args=(u1,params['degradation'], inter))
    p = {}
    for idx, i in enumerate(range(0,24*7+2,1)):
        p[i] = prediction[idx][0]

    res = []
    B_num = 0
    B_den = 0
    for k  in data.keys():
        B_num+=(p[k]*data[k])
        B_den+= p[k]**2
    B = B_num/B_den



    for k  in data.keys():
        res.append((B*p[k]-data[k])**2)
    
    return B
In [35]:
t = numpy.array(t)

pELF3 = read_excel('data/ELF3ELF4LUXProtein.xlsx', sheetname='pELF4')

w_data['cE4'] ={}
for idx,i in enumerate(pELF3.time):
     w_data['cE4'][int(pELF3.iloc[idx,0])+24*5]  =  pELF3.iloc[idx,1]

temp_time, temp = interpol_data(network_data, 'col_0_LL', 'log_cE4_m')
b = tile(temp[0:12],5)
b = concatenate((b,temp[12:]))
p7inter = CubicSpline(range(0,len(b)*2,2),b)
plot(range(0,len(b)*2,2),exp(b))
plot(linspace(2,24*7,200),exp(p7inter(linspace(2,24*7,200))))

figure(dpi=600)
xticks(fontsize=16)
xticks(range(0,24*3,12),fontsize=16)
param = Parameters()
param.add('y',value=1.0173e+05, min=0.0, max=1e6, vary=True)
param.add('degradation',value=0.3848250, min=0.0, max=3.0, vary=True)
out = minimize(residuals,param,args=(w_data['cE4'],p7inter))
report_fit(out.params)
res['ELF4']=odeint(simple_model,array([out.params['y'].value]),t,args=(u1,out.params['degradation'].value,p7inter))
B = sf(out.params, w_data['cE4'],p7inter)
plot(t-24*5,res['ELF4'], 'red', label='pELF4', lw=3)
plot(pELF3.time-24, pELF3.pELF4/B, 'ro-', label='ELF4')
xlabel('Time in constant light (h)', fontsize=22)
xlim(-24,24*2)
#ylim(0,1.4e5)
yticks(fontsize=16)
ylabel('#monomers/cell', fontsize=22)

pro_data['cE4'] = {}
for idx, i in enumerate(res['ELF4']):
    if idx >= 24*4:
        pro_data['cE4'][t[idx]] =  (res['ELF4'][idx][0],0.1)



pELF3 = read_excel('data/ELF3ELF4LUXProtein.xlsx', sheetname='pELF3')
w_data['cE3'] ={}
for idx,i in enumerate(pELF3.time):
     w_data['cE3'][int(pELF3.iloc[idx,0])+24*5]  =  pELF3.iloc[idx,1]

temp_time, temp = interpol_data(network_data, 'col_0_LL', 'log_cE3_m')
b = tile(temp[0:12],5)
b = concatenate((b,temp[12:]))
p7inter = CubicSpline(range(0,len(b)*2,2),b)
#plot(range(0,len(b)*2,2),exp(b))
#plot(linspace(2,24*7,200),exp(p7inter(linspace(2,24*7,200))))


xticks(fontsize=16)
xticks(range(0,24*3,12),fontsize=16)
param = Parameters()
param.add('y',value=1.0173e+05, min=0.0, max=1e6, vary=True)
param.add('degradation',value=0.30848250, min=0.0, max=3.0, vary=True)
out = minimize(residuals,param,args=(w_data['cE3'],p7inter))
report_fit(out.params)
res['ELF3']=odeint(simple_model,array([out.params['y'].value]),t,args=(u1,out.params['degradation'].value,p7inter))
B = sf(out.params, w_data['cE3'],p7inter)
plot(t-24*5,res['ELF3'], 'blue', label='pELF3', lw=3)
plot(pELF3.time-24, pELF3.pELF3/B, 'bo-', label='cE3')
xlabel('Time in constant light (h)', fontsize=22)
xlim(-24,24*2)
#ylim(0,30e3)
yticks(fontsize=16)
ylabel('#monomers/cell', fontsize=22)

pro_data['cE3'] = {}
for idx, i in enumerate(res['ELF3']):
    if idx >= 24*4:
        pro_data['cE3'][t[idx]] =  (res['ELF3'][idx][0],0.1)


pELF3 = read_excel('data/ELF3ELF4LUXProtein.xlsx', sheetname='pLUX')

w_data['cLUX'] ={}
for idx,i in enumerate(pELF3.time):
     w_data['cLUX'][int(pELF3.iloc[idx,0])+24*5]  =  pELF3.iloc[idx,1]

temp_time, temp = interpol_data(network_data, 'col_0_LL', 'log_cLUX_m')
b = tile(temp[0:12],5)
b = concatenate((b,temp[12:]))
p7inter = CubicSpline(range(0,len(b)*2,2),b)
#plot(range(0,len(b)*2,2),exp(b))
#plot(linspace(2,24*7,200),exp(p7inter(linspace(2,24*7,200))))

xticks(fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
xticks(range(0,24*3,12),fontsize=16)
param = Parameters()
param.add('y',value=1.0173e+05, min=0.0, max=1e6, vary=True)
param.add('degradation',value=0.03848250, min=0.0, max=3.0, vary=True)
out = minimize(residuals,param,args=(w_data['cLUX'],p7inter))
report_fit(out.params)
res['LUX']=odeint(simple_model,array([out.params['y'].value]),t,args=(u1,out.params['degradation'].value,p7inter))
B = sf(out.params, w_data['cLUX'],p7inter)

pro_data['cLUX'] = {}
for idx, i in enumerate(res['LUX']):
    if idx >= 24*4:
        print i
        pro_data['cLUX'][t[idx]] =  (res['LUX'][idx][0],0.1)


plot(t-24*5,res['LUX'], 'black', label='pLUX', lw=3)
plot(pELF3.time-24, pELF3.pLUX/B, 'ko-', label='LUX')
xlabel('Time in constant light (h)', fontsize=22)
xlim(-24,24*2)
ylim(1,1e6)
yticks(fontsize=16)
ylabel('#monomers/cell', fontsize=22)
legend(loc='upper right', fontsize=16)
yscale('log')
#savefig('EC_protein_prediction.svg',format='svg',dpi=600,transparent=True)

#savefig('figures/EC_fitting.png', format='png',dpi=300)

[[Variables]]
    y:            101729.934 +/- 18002.0879 (17.70%) (init = 101730)
    degradation:  0.38482477 +/- 0.04124358 (10.72%) (init = 0.384825)
[[Correlations]] (unreported correlations are < 0.100)
    C(y, degradation) = -0.929
[[Variables]]
    y:            107805.104 +/- 6407.37785 (5.94%) (init = 101730)
    degradation:  0.07779919 +/- 0.01510069 (19.41%) (init = 0.3084825)
[[Correlations]] (unreported correlations are < 0.100)
    C(y, degradation) = -0.922
[[Variables]]
    y:            103360.576 +/- 1330503.32 (1287.24%) (init = 101730)
    degradation:  0.04196282 +/- 0.01075398 (25.63%) (init = 0.0384825)
[[Correlations]] (unreported correlations are < 0.100)
    C(y, degradation) =  0.940
[70649.94604933]
[68403.92170012]
[66178.52141009]
[63977.54112552]
[61829.86885218]
[59975.76298822]
[59001.80089699]
[59936.15075344]
[63953.60521998]
[71752.10111323]
[82309.24772231]
[92289.83904187]
[98192.8756314]
[99230.63653349]
[97323.56537845]
[94321.21029063]
[90906.1685419]
[87314.08625754]
[83769.00770673]
[80380.76612246]
[77337.15153228]
[75153.01740683]
[73809.42111158]
[72273.03533385]
[70135.08181061]
[67742.24195775]
[65356.54687234]
[63097.75811601]
[61034.2622367]
[59179.61413707]
[57675.64761174]
[57535.76597543]
[62543.18526579]
[77558.65212044]
[96939.99927261]
[108419.59499932]
[111433.32975803]
[111452.84587023]
[111497.31868486]
[112969.39739504]
[115689.99863066]
[116841.44072113]
[115511.65084306]
[113409.02262869]
[111722.81337717]
[110464.07473148]
[108865.31122066]
[106364.46947618]
[103272.88519688]
[100351.24309622]
[98101.10694738]
[95881.41515551]
[93195.03433159]
[90735.27061503]
[89379.74504352]
[88749.8101671]
[88048.14539937]
[88928.40985501]
[95363.28395825]
[106528.23756302]
[114980.38845825]
[119369.24048392]
[121246.30176596]
[121469.80178563]
[120829.81396236]
[120406.33129398]
[120974.40300529]
[122020.29904494]
[122351.32696723]
[121455.01765552]
[119481.19470812]
[116835.91272031]

Here is were use the new models for rescaling¶

In [36]:
import tellurium as te
from SloppyCell.ReactionNetworks import Utility
In [ ]:
#Berify that the rescaling does not change the model dynamics
In [37]:
%pylab inline
U2017 = te.loada('Models/U2020.4.txt')
U2017res = U2017.simulate(0,24*15,300)
plot(U2017res['time'],U2017res['[cL]']/max(U2017res['[cL]']))
U2017 = te.loada('Models/U2020.4.txt')
U2017.setValue('gl',5)
U2017.setValue('cL',U2017.getValue('cL')*U2017.getValue('gl'))
U2017res = U2017.simulate(0,24*15,300)
#plot(res['time'],res['[cL]'])

plot(U2017res['time'],U2017res['[cL]']/max(U2017res['[cL]']))

Populating the interactive namespace from numpy and matplotlib

/usr/local/lib/python2.7/dist-packages/IPython/core/magics/pylab.py:161: UserWarning:

pylab import has clobbered these variables: ['plotting', 'array', 'datetime', 'test', 'unique']
`%matplotlib` prevents importing * from pylab and numpy
Out[37]:
[<matplotlib.lines.Line2D at 0x40a079c4d0>]
In [38]:
U2017 = te.loada('Models/U2020.4.txt')
U2017res = U2017.simulate(0,24*15,300)
plot(U2017res['time'],U2017res['[cP7]']/max(U2017res['[cP7]']))
U2017 = te.loada('Models/U2020.4.txt')
U2017.setValue('gp7',2)
U2017.setValue('cP9',U2017.getValue('cP9')*U2017.getValue('gp7'))
U2017.setValue('cP7',U2017.getValue('cP7')*U2017.getValue('gp7'))
U2017.setValue('cP5',U2017.getValue('cP5')*U2017.getValue('gp7'))
U2017.setValue('cT',U2017.getValue('cT')*U2017.getValue('gp7'))
res = U2017.simulate(0,24*10,300)
plot(U2017res['time'],U2017res['[cP7]']/max(U2017res['[cP7]']))
Out[38]:
[<matplotlib.lines.Line2D at 0x40b16888d0>]
In [39]:
%pylab inline
U2017 = te.loada('Models/U2020.4.txt')
#
U2017.reset() #= te.loada('Models/U2017.txt')


U2017.setValue('ge3',1)
range(0,24*15,1)
res = U2017.simulate(0,24*15,24*15)
#plot(res['time'],res['[cLUX]']/max(res['[cLUX]']))
plot(res['time'],res['[cEC]']/max(res['[cEC]']))
U2017 = te.loada('Models/U2020.4.txt')

U2017.setValue('ge3',100)
#xlim(-24,24*5)
xticks(range(0,24*5,24))
U2017.setValue('cE3c',U2017.getValue('cE3c')*U2017.getValue('ge3'))
U2017.setValue('cE3n',U2017.getValue('cE3n')*U2017.getValue('ge3'))
U2017.setValue('cE4',U2017.getValue('cE4')*U2017.getValue('ge3'))
U2017.setValue('cEGc',U2017.getValue('cEGc')*U2017.getValue('ge3'))
U2017.setValue('cEGn',U2017.getValue('cEGn')*U2017.getValue('ge3'))
U2017.setValue('cLUX',U2017.getValue('cLUX')*U2017.getValue('ge3'))
U2017.setValue('cGc',U2017.getValue('cGc')*U2017.getValue('ge3'))
U2017.setValue('cGn',U2017.getValue('cGn')*U2017.getValue('ge3'))
U2017.setValue('cEC',U2017.getValue('cEC')*U2017.getValue('ge3'))
U2017.setValue('cZTL',U2017.getValue('cZTL')*U2017.getValue('ge3'))
U2017.setValue('cZG',U2017.getValue('cZG')*U2017.getValue('ge3'))
U2017.setValue('cE34',U2017.getValue('cE34')*U2017.getValue('ge3'))
res = U2017.simulate(0,24*15,300)
#plot(res['time'],res['[cLUX]']/max(res['[cLUX]']))
plot(res['time'],res['[cEC]']/max(res['[cEC]']))
#plot(res['time'],res['[cLUX]'])

Populating the interactive namespace from numpy and matplotlib
Out[39]:
[<matplotlib.lines.Line2D at 0x40a0a3be90>]
In [40]:
m = pro_data['cL'].values()[0][0] 
for i in pro_data['cL'].values():
    if i[0] > m:
        m = i[0]
m
Out[40]:
106759.4287181982
In [41]:
U2017.resetAll()
def residuals(params, U2017, id,gene, pro_data):
    U2017.resetAll()
    U2017.setValue(id, params[id].value)
    U2017.setValue(gene,U2017.getValue(gene)*params[id].value)
    res = U2017.simulate(0,24*15,24*15)
    
    max(res['['+gene+']'])
    
    residuals = []
    pred = {}
    for idx,i in enumerate(res['['+gene+']']):
        pred[around(res['time'][idx])-24*5] = res['['+gene+']'][idx]
    for i in pro_data[gene].keys():
            residuals.append(((pro_data[gene][i][0])-(pred[int(i)]))**2)
    return residuals

def residuals_max(params, U2017, id,gene, pro_data):
    U2017.resetAll()
    U2017.setValue(id, params[id].value)
    U2017.setValue(gene,U2017.getValue(gene)*params[id].value)
    res = U2017.simulate(0,24*15,24*15)
    
    max(res['['+gene+']'])
    
    residuals = []
    pred = {}
    for idx,i in enumerate(res['['+gene+']']):
        pred[around(res['time'][idx])-24*5] = res['['+gene+']'][idx]
    
    m1 = pro_data[gene].values()[0][0] 
    for i in pro_data[gene].values():
        if i[0] > m1:
            m1 = i[0]
    
    for i in pro_data[gene].keys():
            residuals.append(((pro_data[gene][i][0])-max(res['['+gene+']']))**2)
    return (m1-max(res['['+gene+']']))**2
In [42]:
p = {}
params = Parameters()
params.add('gl',value=1e5, min=1e4, max=1e6, vary=True)
out = minimize(residuals_max,params,args=(U2017,'gl','cL',pro_data), method='leastsq')
report_fit(out.params)
p['gl'] = out.params['gl'].value

[[Variables]]
    gl:  99593.2079 +/- 197.600037 (0.20%) (init = 100000)
In [44]:
figure(dpi=600)
lh = []
r = linspace(1,0.2e6,100)
for i in r:
    params['gl'].value=i
    lh.append(sum(residuals_max(params, U2017, 'gl','cL', pro_data)))

yticks(fontsize=16)
xticks(fontsize=16)
plot(r,lh)
ylabel('Cost',fontsize=22)
xlabel('scale factor value',fontsize=22)
title('Cost landscape', fontsize=22)
tight_layout()
xscale('log')
savefig('cost_landscape_cL_max_U2020_4.svg', format='svg',dpi=600, transparent=True)
In [45]:
params = Parameters()
params.add('gp7',value=2e4, min=1e4, max=1e6, vary=True)
out = minimize(residuals,params,args=(U2017,'gp7','cP7',pro_data), method='leastsq')
report_fit(out.params)
p['gp7'] = out.params['gp7'].value

[[Variables]]
    gp7:  59772.0066 +/- 34.5369483 (0.06%) (init = 60000)
In [46]:
figure(dpi=600)
lh = []
r = linspace(1,0.8e5,100)
for i in r:
    params['gp7'].value=i
    lh.append(sum(residuals(params, U2017, 'gp7','cP7', pro_data)))

plot(r,lh)

yticks(fontsize=16)
xticks(fontsize=16)
plot(r,lh)
ylabel('Cost',fontsize=22)
xlabel('scale factor value',fontsize=22)
title('Cost landscape', fontsize=22)
tight_layout()
xscale('log')
savefig('cost_landscape_P7_max_U2020_4.svg', format='svg',dpi=600, transparent=True)
In [47]:
params = Parameters()
params.add('gp7',value=2e4, min=1e4, max=1e6, vary=True)
out = minimize(residuals,params,args=(U2017,'gp7','cP7',pro_data), method='leastsq')
report_fit(out.params)
p['gp7'] = out.params['gp7'].value

[[Variables]]
    gp7:  19558.4309 +/- 82.0318414 (0.42%) (init = 20000)
In [48]:
def residuals(params, U2017, id,gene, pro_data):
    U2017.resetAll()
    #for i in U2017p.keys():
    #    U2017.setValue(i,U2017p.getByKey(i))
    U2017.setValue(id, params[id].value)
    U2017.setValue(gene,U2017.getValue(gene)*params[id].value)
    U2017.setValue('cE3c',U2017.getValue('cE3c')*params[id].value)
    U2017.setValue('cE3n',U2017.getValue('cE3n')*params[id].value)
    U2017.setValue('cE4',U2017.getValue('cE4')*params[id].value)
    U2017.setValue('cEGc',U2017.getValue('cEGc')*params[id].value)
    U2017.setValue('cEGn',U2017.getValue('cEGn')*params[id].value)
    #U2017.setValue('cLUX',U2017.getValue('cLUX')*params[id].value)
    U2017.setValue('cGc',U2017.getValue('cGc')*params[id].value)
    U2017.setValue('cGn',U2017.getValue('cGn')*params[id].value)
    U2017.setValue('cEC',U2017.getValue('cEC')*params[id].value)
    U2017.setValue('cZTL',U2017.getValue('cZTL')*params[id].value)
    U2017.setValue('cZG',U2017.getValue('cZG')*params[id].value)
    U2017.setValue('cE34',U2017.getValue('cE34')*params[id].value)
    res = U2017.simulate(0,24*15,24*15)
    residuals = []
    pred = {}
    for idx,i in enumerate(res['['+gene+']']):
        pred[around(res['time'][idx])-24*5] = res['['+gene+']'][idx]
    for i in pro_data[gene].keys():
            residuals.append(((pro_data[gene][i][0])-(pred[int(i)]))**2)
    return residuals
In [55]:
params = Parameters()
params.add('ge3',value=1e4, min=1e2, max=1e6, vary=True)
out = minimize(residuals,params,args=(U2017,'ge3','cLUX',pro_data), method='leastsq')
report_fit(out.params)
params = out.params
p['ge3'] = out.params['ge3'].value

[[Variables]]
    ge3:  10277.7345 +/- 39.4261097 (0.38%) (init = 10000)
In [56]:
figure(dpi=600)
lh = []
r = linspace(1e3,5e4,100)
for i in r:
    params['ge3'].value=i
    lh.append(sum(residuals(params, U2017, 'ge3','cLUX', pro_data)))

plot(r,lh)

yticks(fontsize=16)
xticks(fontsize=16)
plot(r,lh)
ylabel('Cost',fontsize=22)
xlabel('scale factor value',fontsize=22)
title('Cost landscape', fontsize=22)
tight_layout()
xscale('log')
savefig('cost_landscape_E3_max_U2020_4.svg', format='svg',dpi=600, transparent=True)
In [58]:
p
Out[58]:
{'ge3': 10277.734457469609, 'gl': 99593.20789266606, 'gp7': 19558.430910820854}
In [59]:
U2017_scaling_factors = DataFrame.from_dict(p, orient="index")
U2017_scaling_factors.to_excel("U2020_4_scaling_factors.xls")
In [60]:
p_loaded = read_excel("U2020_4_scaling_factors.xls")
In [61]:
sf={}
for idx, scaling_factor in enumerate(p_loaded.iloc[:,0]):
    sf[scaling_factor]= p_loaded.iloc[idx,1]
sf
Out[61]:
{u'ge3': 10277.734457469609,
 u'gl': 99593.20789266606,
 u'gp7': 19558.430910820854}
In [62]:
U2017 = te.loada('Models/U2020.4.txt')

for i in p.keys():
    U2017.setValue(i,sf[i])
#U2017.setValue('m15',0.3)
res = U2017.simulate(0,24*15,3000)
In [75]:
U2017.exportToAntimony("Models/U2020.4.scaled.txt")
U2017.exportToSBML("Models/U2020.4.scaled.sbml")
In [64]:
get_max = []
figure(dpi=600)
plot(res['time']-24*11,(res['[cL]']), color='k')
xlim(-24,24*2)
yticks(fontsize=16)
xticks(range(0,24*2,12), fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(24*2+12,24*3, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
d = []
t = []
for i in pro_data['cL'].keys():
    get_max.append(pro_data['cL'][i][0])
    d.append(pro_data['cL'][i][0])
    t.append(i-24*5)
    
d = array(d)
t = array(t)
t = concatenate((t[range(41,72)],t[range(0,40)]))
d = concatenate((d[range(41,72)],d[range(0,40)]))

plot(t,d,'-', color='orange', lw=3)

ylabel('#monomers/cell', fontsize=22)
xlabel('Time hours in constant light (h)', fontsize=22)
yscale('log')
ylim(1e1,1e6)
tight_layout()

savefig('U2020_4_cL_rescaled_total.svg',format='svg',dpi=600, transparent=True)
print max(get_max[65:])
print min(get_max[65:])

70339.46919568657
4848.913455604003
In [92]:
def plot_dict(data,var,sf=1):
    time = array(sorted(data[var].keys()))
    d = []
    for i in sorted(data[var].keys()):
        d.append(data[var][i])

    plot(time-24*5,array(d)*sf,color='blue')
In [66]:
figure(dpi=600)
plot(res['time']-24*11,(res['[cP9]']), 'k')
d = []
t = []
for i in pro_data['cP9'].keys():
    #plot(i-24*5,pro_data['cP9'][i][0],'o', color='orange')
    d.append(pro_data['cP9'][i][0])
    t.append(i-24*5)

d = array(d)
t = array(t)
t = concatenate((t[range(41,72)],t[range(0,40)]))
d = concatenate((d[range(41,72)],d[range(0,40)]))    

plot(t,d,'-', color='orange', lw=3)   
#plot_dict(w_data,'cP9',sf=10000)
yticks(fontsize=16)
xlim(-24,24*2)
ylim(1,1e6)
xticks(range(0,24*2,12),fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
xlabel('Time in constant light (h)', fontsize=22)
ylabel('#monomers/cell', fontsize=22)
yscale('log')
tight_layout()

savefig('U2020_4_cP9_rescaled_total.svg',format='svg',dpi=600, transparent=True)
In [67]:
figure(dpi=600)
plot(res['time']-24*11,(res['[cP7]']), 'k')
d = []
t = []
for i in pro_data['cP7'].keys():
    #plot(i-24*5,pro_data['cP9'][i][0],'o', color='orange')
    d.append(pro_data['cP7'][i][0])
    t.append(i-24*5)

d = array(d)
t = array(t)
t = concatenate((t[range(41,72)],t[range(0,40)]))
d = concatenate((d[range(41,72)],d[range(0,40)]))    

plot(t,d,'-', color='orange', lw=3)   
#plot_dict(w_data,'cP9',sf=10000)
yticks(fontsize=16)
xlim(-24,24*2)
ylim(1e2,1e5)
xticks(range(0,24*2,12),fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
xlabel('Time in constant light (h)', fontsize=22)
ylabel('#monomers/cell', fontsize=22)
yscale('log')
tight_layout()
savefig('U2020_4_cP7_rescaled_total.svg',format='svg',dpi=600, transparent=True)
In [68]:
figure(dpi=600)
plot(res['time']-24*11,(res['[cP5]']), 'k')
d = []
t = []
for i in pro_data['cP5'].keys():
    #plot(i-24*5,pro_data['cP9'][i][0],'o', color='orange')
    d.append(pro_data['cP5'][i][0])
    t.append(i-24*5)

d = array(d)
t = array(t)
t = concatenate((t[range(41,72)],t[range(0,40)]))
d = concatenate((d[range(41,72)],d[range(0,40)]))    

plot(t,d,'-', color='orange', lw=3)   
#plot_dict(w_data,'cP9',sf=10000)
yticks(fontsize=16)
xlim(-24,24*2)
ylim(1e2,1e5)
xticks(range(0,24*2,12),fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
xlabel('Time in constant light (h)', fontsize=22)
ylabel('#monomers/cell', fontsize=22)
yscale('log')
tight_layout()
savefig('U2020_4_cP5_rescaled_total.svg',format='svg',dpi=600, transparent=True)
In [69]:
figure(dpi=600)
plot(res['time']-24*11,(res['[cT]']), 'k')
d = []
t = []
for i in pro_data['cT'].keys():
    #plot(i-24*5,pro_data['cP9'][i][0],'o', color='orange')
    d.append(pro_data['cT'][i][0])
    t.append(i-24*5)

d = array(d)
t = array(t)
t = concatenate((t[range(41,72)],t[range(0,40)]))
d = concatenate((d[range(41,72)],d[range(0,40)]))    

plot(t,d,'-', color='orange', lw=3)   
#plot_dict(w_data,'cP9',sf=10000)
yticks(fontsize=16)
xlim(-24,24*2)
ylim(1e2,1e5)
xticks(range(0,24*2,12),fontsize=16)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
xlabel('Time in constant light (h)', fontsize=22)
ylabel('#monomers/cell', fontsize=22)
yscale('log')
tight_layout()
savefig('U2020_4_cT_rescaled_total.svg',format='svg',dpi=600, transparent=True)
In [98]:
## print log(2)/U2017.m15
print log(2)/U2017.m23
print log(2)/U2017.m17
print log(2)/U2017.m24

0.4681601576364287
1.1426575868077768
9.040130030680263
In [208]:
temp_pro_data = []
temp_pro_data_time = []
for i in sorted(pro_data['cT'].keys()):
    temp_pro_data.append(pro_data['cT'][i][0])
    temp_pro_data_time.append(i)
In [209]:
#plot(temp_pro_data_time[:20],temp_pro_data[:20]/max(temp_pro_data[:20]),'o')
plot(temp_pro_data_time,temp_pro_data,'o')
print min(temp_pro_data[20:])
print max(temp_pro_data[20:])

452.94743129192386
29981.58067324995
In [151]:
figure(figsize=(10,6))
yticks(fontsize=16)
xticks(range(0,24*2,12),fontsize=16)
plot(res['time']-24*11,(res['[cP7]'])/max(res['[cP7]']),'k', label='cP7')
plot(res['time']-24*11,(res['[cL_m]']*150)/max(res['[cL_m]']*150),'b', label='cL_m')
rep = U2017.gp7**2*U2017.g1**2/(U2017.gp7**2*U2017.g1**2+res['[cP7]']**2)
plot(res['time']-24*11,1-rep, 'r', label='cP7 repression')
plot(array(temp_pro_data_time)-24*5,temp_pro_data/max(temp_pro_data),'o', color='orange', label='PRR7 SPD')
#plot_dict(w_data,'cP7',sf=50000)
title('cP7', fontsize=22)
xlim(-24,24*2)
ylim(0,1.7)
axvspan(-12,0, color='gray',alpha=0.7)
ylabel('Normalised to max', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
axvspan(24+12,24*2, color='gray',alpha=0.1)
axvspan(12,24, color='gray',alpha=0.1)
legend(loc='upper right', fontsize=16)

temp_time =[]
temp_rna =[]
for i in sorted(network_data['col_0_LL']['log_cC_m'].keys()):
    temp_time.append(i)
    temp_rna.append(exp(network_data['col_0_LL']['log_cC_m'][i][0]))

plot(array(temp_time[:])-24*10,temp_rna[:]/max(temp_rna[:]),'bo-')

savefig('cL_rise.png', format='png',dpi=300)
In [152]:
temp_time =[]
temp_rna =[]
for i in sorted(network_data['col_0_LL']['log_cC_m'].keys()):
    temp_time.append(i)
    temp_rna.append(exp(network_data['col_0_LL']['log_cC_m'][i][0]))

plot(temp_time[1:],temp_rna[1:],'bo-')
Out[152]:
[<matplotlib.lines.Line2D at 0x12bb3a610>]
In [153]:
TOPLESS_expression = array([337.836,266.11,225.958,399.275,470.614,438.132,
                            320.583,342.415,279.855,428.996,411.956,324.224])
plot(range(0,24*2,4),TOPLESS_expression/max(TOPLESS_expression),'bo-')
Out[153]:
[<matplotlib.lines.Line2D at 0x129fd1c50>]
In [154]:
figure(figsize=(10,6))
yticks(fontsize=16)
xticks(range(0,24*2,12),fontsize=16)
plot(res['time']-24*11,(res['[cP7]'])/max(res['[cP7]']),'k', label='cP7')
plot(res['time']-24*11,(res['[cL_m]']*150)/max(res['[cL_m]']*150),'b', label='cL_m')
plot(range(0,24*2,4),TOPLESS_expression/max(TOPLESS_expression),'go-', label='TPL Diurnal Database')
plot(array(temp_pro_data_time)-24*5,temp_pro_data/max(temp_pro_data),'o', color='orange', label='PRR7 SPD')
#plot_dict(w_data,'cP7',sf=50000)
title('cP7', fontsize=22)
xlim(-24,24*2)
ylim(0,1.7)
axvspan(-12,0, color='gray',alpha=0.7)
ylabel('Normalised to max', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
axvspan(24+12,24*2, color='gray',alpha=0.1)
axvspan(12,24, color='gray',alpha=0.1)
legend(loc='upper right', fontsize=16)
savefig('topless_progile.png', format='png',dpi=300)
In [155]:
U2017sf
Out[155]:
KeyedList([('gml', 9696.3284783153849), ('gm9', 3645.1361509496446), ('gm7', 1854.8827904077939), ('gm5', 530.27711530965303), ('gmt', 4881.4791683880294), ('gm4', 4919.2311439361465), ('gm3', 2847.6321058166104), ('gmg', 4987.7906816742825)])
In [211]:
U2017k = te.loada('Models/U2017kinase.txt')

U2017p = Utility.load('Models/params_4.bp')
U2017sf = Utility.load('Models/sfU2017.bp')
for i in U2017p.keys():
    U2017k.setValue(i,U2017p.getByKey(i))
for i in U2017sf.keys():
    U2017k.setValue(i,U2017sf.getByKey(i))
for i in p.keys():
    U2017k.setValue(i,p[i])
U2017k.setValue('k_strength',0)
reskbefore = U2017k.simulate(0,24*15,3000)
print log(2)/U2017k.m15
print log(2)/U2017k.m23
U2017k.setValue('k_strength',0.15)
U2017k.setValue('m15',0.3)
U2017k.setValue('m23',0.01)
print log(2)/0.3
print log(2)/0.01
resk = U2017k.simulate(0,24*15,3000)

1.004126352466796
0.39976022536822303
2.3104906018664844
69.31471805599453
In [216]:
figure(figsize=(10,6))
yticks(fontsize=16)
xticks(range(0,24*2,12),fontsize=16)
plot(res['time']-24*11,(resk['[cP7]']),'k', label='cP7')
plot(res['time']-24*11,(resk['[cL_m]']*150),'b', label='cL_m*150')
plot(reskbefore['time']-24*11,(reskbefore['[cL_m]']*150),'b--', label='cL_m*150 no K')
plot(res['time']-24*11,(resk['[cK]']*150),'r', label='HKinase')
#plot(range(0,24*2,4),TOPLESS_expression,'go-', label='TPL Diurnal Database')
#plot(array(temp_pro_data_time)-24*5,temp_pro_data,'o', color='orange', label='PRR7 SPD')
#plot_dict(w_data,'cP7',sf=50000)
#title('cP7', fontsize=22)
xlim(-24,24*2)
#ylim(0,1.7)
axvspan(-12,0, color='gray',alpha=0.7)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
axvspan(12+24*2,24*3, color='gray',alpha=0.3)
temp_time =[]
temp_rna =[]
for i in sorted(network_data['col_0_LL']['log_cC_m'].keys()):
    temp_time.append(i)
    temp_rna.append(exp(network_data['col_0_LL']['log_cC_m'][i][0]))

plot(array(temp_time[:])-24*11,array(temp_rna[:])*1.3e2,'bo-', label='CCA1 mRNA TiMet')
legend(loc='upper right', fontsize=16)
yscale('log')
#savefig('hkinase_profile.png', format='png',dpi=300)
In [71]:
figure(dpi=600)
d = []
t = []
for i in pro_data['cE3'].keys():
    #plot(i-24*5,pro_data['cP9'][i][0],'o', color='orange')
    d.append(pro_data['cE3'][i][0])
    t.append(i-24*5)

d = array(d)
t = array(t)
t = concatenate((t[range(41,72)],t[range(0,40)]))
d = concatenate((d[range(41,72)],d[range(0,40)]))   
plot(t,d,'-', color='orange', lw=3) 

xticks(range(0,24*2,12),fontsize=16)
plot(res['time']-24*11,(res['[cE3c]']), 'k')
#for i in pro_data['cE3'].keys():
#    plot(i-24*5,pro_data['cE3'][i][0],'o', color='orange')
#plot_dict(w_data,'cE3',sf=30000)
yticks(fontsize=16)
xlim(-24,24*2)
ylim(1e1,1e5)
yscale('log')
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
tight_layout()
savefig('U2020_4_cE3_rescaled.svg',format='svg',dpi=600,transparent=True)
In [72]:
figure(dpi=600)
d = []
t = []
for i in pro_data['cE4'].keys():
    #plot(i-24*5,pro_data['cP9'][i][0],'o', color='orange')
    d.append(pro_data['cE4'][i][0])
    t.append(i-24*5)

d = array(d)
t = array(t)
t = concatenate((t[range(41,72)],t[range(0,40)]))
d = concatenate((d[range(41,72)],d[range(0,40)]))   
plot(t,d,'-', color='orange', lw=3) 

xticks(range(0,24*2,12),fontsize=16)
plot(res['time']-24*11,(res['[cE4]']), 'k')
#for i in pro_data['cE3'].keys():
#    plot(i-24*5,pro_data['cE3'][i][0],'o', color='orange')
#plot_dict(w_data,'cE3',sf=30000)
yticks(fontsize=16)
xlim(-24,24*2)
ylim(1e2,1e6)
yscale('log')
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
tight_layout()
savefig('U2020_4_cE4_rescaled.svg',format='svg',dpi=600,transparent=True)
In [74]:
figure(dpi=600)
d = []
t = []
for i in pro_data['cLUX'].keys():
    #plot(i-24*5,pro_data['cP9'][i][0],'o', color='orange')
    d.append(pro_data['cLUX'][i][0])
    t.append(i-24*5)

d = array(d)
t = array(t)
t = concatenate((t[range(41,72)],t[range(0,40)]))
d = concatenate((d[range(41,72)],d[range(0,40)]))   
plot(t,d,'-', color='orange', lw=3) 

xticks(range(0,24*2,12),fontsize=16)
plot(res['time']-24*11,(res['[cLUX]']), 'k')
#for i in pro_data['cE3'].keys():
#    plot(i-24*5,pro_data['cE3'][i][0],'o', color='orange')
#plot_dict(w_data,'cE3',sf=30000)
yticks(fontsize=16)
xlim(-24,24*2)
ylim(1e3,1e6)
yscale('log')
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
tight_layout()
savefig('U2020_4_cLUX_rescaled.svg',format='svg',dpi=600,transparent=True)

Translation rates

In [167]:
print max(get_max[20:])
print min(get_max[20:])

128663.174159
62252.9310291
In [240]:
pbm_affinities = read_excel('data/LUX_affinities_PBM_U2017.xlsx')
pbm_affinities
Out[240]:
Model Gene promoter Kd value
0 PBM LUX 0.260289 0.260289
1 PBM PRR9 0.276830 0.276830
2 PBM CCA1 0.393578 0.393578
3 PBM PRR5 1.977416 1.977416
4 PBM TOC1 0.685716 0.685716
5 U2017 PRR9 9.250000 0.010300
6 U2017 TOC1 10.200000 0.011500
7 U2017 LUX 8.850000 0.009900
8 U2017 GI 3.580000 0.004010
In [241]:
plt.style.use('seaborn-white')
figure(figsize=(10,6))
sns.set_palette("Paired")
pbm_vs_U2017 =  read_excel('data/LUX_affinities_PBM_U2017.xlsx')
sns.barplot(x='Gene promoter', y='Kd',data=pbm_vs_U2017, hue='Model')
legend(loc='upper left', fontsize=16)
ylabel('$k_d$ nM', fontsize=22)
xlabel('Gene promoter',fontsize=22)
ylim(0,12)
#yscale('log')
xticks(fontsize=16)
yticks(fontsize=16)
savefig('EC_Kd_comparision_absolute.png', format='png', dpi=300)


NameErrorTraceback (most recent call last)
<ipython-input-241-38351d91821a> in <module>()
      1 plt.style.use('seaborn-white')
      2 figure(figsize=(10,6))
----> 3 sns.set_palette("Paired")
      4 pbm_vs_U2017 =  read_excel('data/LUX_affinities_PBM_U2017.xlsx')
      5 sns.barplot(x='Gene promoter', y='Kd',data=pbm_vs_U2017, hue='Model')

NameError: name 'sns' is not defined
<Figure size 720x432 with 0 Axes>
In [242]:
U2017k = te.loada('Models/U2017KBOA.txt')

U2017p = Utility.load('Models/params_4.bp')
U2017sf = Utility.load('Models/sfU2017.bp')
for i in U2017p.keys():
    U2017k.setValue(i,U2017p.getByKey(i))
for i in U2017sf.keys():
    U2017k.setValue(i,U2017sf.getByKey(i))
for i in p.keys():
    U2017k.setValue(i,p[i])


U2017k.setValue('k_strength',0.1)
U2017k.setValue('pBOA',0.0)
U2017k.setValue('pBOA',20)

U2017k.setValue('pBOA2',0.5)
U2017k.setValue('pBOA2',200)
#U2017k.setValue('p27',0)
U2017k.setValue('m15',0.3)
U2017k.setValue('m23',0.3)
#U2017k.setValue('aff',0.1)
U2017k.setValue('aff',0.01)
U2017k.setValue('fass_ev',200)
#U2017k.setValue('fass_ev',1)
print log(2)/0.3
print log(2)/0.01
resk = U2017k.simulate(0,24*15,3000)

figure(figsize=(10,6))
yticks(fontsize=16)
xticks(range(0,24*4,12),fontsize=16)
plot(res['time']-24*11,(resk['[cBOA]']),'g', label='cBOA')
#plot(res['time']-24*11,(resk['[cEC]']),'r-', label='cEC')
plot(res['time']-24*11,(resk['[cL_m]']),'b-', label='BOA-OX cL_m')
xlim(-24,24*10)
#ylim(0,1.7)
axvspan(-12,0, color='gray',alpha=0.7)
ylabel('log(#molecules/cell)', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
axvspan(24+12,24*2, color='gray',alpha=0.1)
axvspan(12,24, color='gray',alpha=0.1)
axvspan(12+24*2,24*3, color='gray',alpha=0.1)
yscale('log')


title('B', fontsize=22)
legend(loc='upper right', fontsize=16)
#savefig('BOA_intro.png',format='png',dpi=300)

2.3104906018664844
69.31471805599453
Out[242]:
<matplotlib.legend.Legend at 0x7f5e359a8910>
In [274]:
amin_sec = 3.
rib_tranc = 10.
PRR7_length = 468.
PRR7_trans_rate = 1/(PRR7_length/amin_sec/10)*3600
plot(res['time']-24*11,res['[cP9_m]']*PRR7_trans_rate,'k', label='Theoretical Max trans')
plot(res['time']-24*11,U2017.getValue('gp7')*U2017.getValue('p8')*res['[cP9_m]']/U2017.getValue('gm9'), 'b--',label='Rescaled')
xticks(range(0,24*2,12), fontsize=22)
yticks(fontsize=22)
xlim(-24,24*2)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
legend(loc='lower right', fontsize=16)
ylabel('#molecules/h', fontsize=22)
yscale('log')
In [270]:
figure(figsize=(12,10), dpi=600)

subplot(221)
amin_sec = 3.
rib_tranc = 10.
PRR7_length = 468.
PRR7_trans_rate = 1/(PRR7_length/amin_sec/10)*3600
plot(res['time']-24*11,res['[cP9_m]']*PRR7_trans_rate,'k', label='Theoretical Max trans')
plot(res['time']-24*11,U2017.getValue('gp7')*U2017.getValue('p8')*res['[cP9_m]']/U2017.getValue('gm9'), 'b--',label='Rescaled')
yscale('log')
xticks(range(0,24*2,12), fontsize=22)
yticks(fontsize=22)
xlim(-24,24*2)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
legend(loc='lower right', fontsize=16)
ylabel('#molecules/h', fontsize=22)

title('cP9', fontsize=22)

subplot(222)
amin_sec = 3.
rib_tranc = 10.
PRR7_length = 618.
PRR7_trans_rate = 1/(PRR7_length/amin_sec/10)*3600
plot(res['time']-24*11,res['[cP7_m]']*PRR7_trans_rate,'k', label='Theoretical max trans.')
plot(res['time']-24*11,U2017.getValue('gp7')*U2017.getValue('p9')*res['[cP7_m]']/U2017.getValue('gm7'), 'b--',label='Rescaled')
yscale('log')
xticks(range(0,24*2,12), fontsize=22)
yticks(fontsize=22)
xlim(-24,24*2)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)


title('cP7', fontsize=22)

subplot(223)
amin_sec = 3.
rib_tranc = 10.
PRR7_length = 558.
PRR7_trans_rate = 1/(PRR7_length/amin_sec/10)*3600
plot(res['time']-24*11,res['[cP5_m]']*PRR7_trans_rate,'k', label='Theoretical Max trans')
plot(res['time']-24*11,U2017.getValue('gp7')*U2017.getValue('p10')*res['[cP5_m]']/U2017.getValue('gm5'),'b--', label='U2017 trans')
yscale('log')
xticks(range(0,24*2,12), fontsize=22)
yticks(fontsize=22)
xlim(-24,24*2)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)

ylabel('#molecules/h', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
title('cP5', fontsize=22)

subplot(224)
amin_sec = 3.
rib_tranc = 10.
PRR7_length = 618.
PRR7_trans_rate = 1/(PRR7_length/amin_sec/10)*3600
plot(res['time']-24*11,res['[cT_m]']*PRR7_trans_rate,'k', label='Theoretical Max trans')
plot(res['time']-24*11,U2017.getValue('gp7')*U2017.getValue('p4')*res['[cT_m]']/U2017.getValue('gmt'), 'b--',label='Rescaled')
yscale('log')
xticks(range(0,24*2,12), fontsize=22)
yticks(fontsize=22)
xlim(-24,24*2)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)


xlabel('Time in constant light (h)', fontsize=22)
title('cT', fontsize=22)
tight_layout()
savefig('PRRS__translation_rates.svg', format='svg', dpi=600, transparent=True)
In [259]:
from scipy.constants import Avogadro
In [260]:
print U2017.getValue('g9')
mean(array([U2017.getValue('g9'),
            U2017.getValue('g10'),
            U2017.getValue('g5'),
            U2017.getValue('g6'),
            U2017.getValue('g16'),
            U2017.getValue('g15'),
            U2017.getValue('g12')]))

0.267310750786
Out[260]:
0.29790048402767483
In [261]:
print 'g9\t', U2017.getValue('gl')*(U2017.getValue('p1')+U2017.getValue('p2'))/U2017.getValue('gml')

g9	158.249727823
In [262]:
U2017.getValue('m3')*2
Out[262]:
0.40110991504829296
In [263]:
A = 0.1168e-12
B = 0.112533e-12
(U2017.getValue('g9')*U2017.getValue('gl')/1100)/(Avogadro*.112533e-12)
Out[263]:
3.0448506487979457e-10
In [264]:
print 'g9\t', U2017.getValue('g9')*U2017.getValue('gl')/(Avogadro*.112533e-12)*1e9
print 'g10\t', U2017.getValue('g10')*U2017.getValue('gl')/(Avogadro*0.112533e-12)*1e9
print 'g12\t', U2017.getValue('g12')*U2017.getValue('gl')/(Avogadro*0.112533e-12)*1e9
print 'g5\t', U2017.getValue('g5')*U2017.getValue('gl')/(Avogadro*0.112533e-12)*1e9
print 'g6\t', U2017.getValue('g6')*U2017.getValue('gl')/(Avogadro*0.112533e-12)*1e9
print 'g16\t', U2017.getValue('g16')*U2017.getValue('gl')/(Avogadro*0.112533e-12)*1e9
print 'g15\t', U2017.getValue('g15')*U2017.getValue('gl')/(Avogadro*0.112533e-12)*1e9

g9	334.933571368
g10	619.642572909
g12	271.418814406
g5	179.422681504
g6	353.102873517
g16	383.312968504
g15	470.998621027
In [265]:
PBM_affinties = array([0.34,
0.25,
0.57,
0.44,
0.70,
0.40,
0.31,
0.40])



print PBM_affinties

PBM_affinties = array([0.44,
0.31,
0.47,
0.4,
0.25,
0.34])

U2017_cca1_aff  = array([U2017.getValue('g9')*U2017.getValue('gl')/(Avogadro*A)*1e9/1100,
U2017.getValue('g10')*U2017.getValue('gl')/(Avogadro*A)*1e9/1100,
U2017.getValue('g12')*U2017.getValue('gl')/(Avogadro*A)*1e9/1100,
U2017.getValue('g5')*U2017.getValue('gl')/(Avogadro*A)*1e9/1100,
U2017.getValue('g6')*U2017.getValue('gl')/(Avogadro*A)*1e9/1100,
U2017.getValue('g16')*U2017.getValue('gl')/(Avogadro*A)*1e9/1100,
U2017.getValue('g15')*U2017.getValue('gl')/(Avogadro*A)*1e9/1000])

U2017_cca1_aff  = U2017_cca1_aff
print U2017_cca1_aff*1100
print U2017_cca1_aff
print PBM_affinties

[0.34 0.25 0.57 0.44 0.7  0.4  0.31 0.4 ]
[322.6975992  597.00545939 261.50319727 172.86791625 340.2031307
 369.30957436 499.17103255]
[0.29336145 0.54273224 0.23773018 0.15715265 0.30927557 0.33573598
 0.45379185]
[0.44 0.31 0.47 0.4  0.25 0.34]
In [266]:
import seaborn as sns


ImportErrorTraceback (most recent call last)
<ipython-input-266-a84c0541e888> in <module>()
----> 1 import seaborn as sns

ImportError: No module named seaborn
In [267]:
pbm_vs_U2017 =  read_excel('data/PBM_vs_U20172.xlsx')
pbm_vs_U2017
Out[267]:
Gene promoter Model value ratio
0 PRR9 U2017.2 0.440000 0.909091
1 PRR7 U2017.2 0.310000 1.290323
2 PRR5 U2017.2 0.470000 0.851064
3 TOC1 U2017.2 0.400000 1.000000
4 LUX U2017.2 0.340000 1.176471
5 ELF4 U2017.2 0.250000 1.600000
6 GI U2017.2 0.570000 0.701754
7 PRR9 EMA selected 0.293360 0.535696
8 PRR7 EMA selected 0.542730 0.289558
9 PRR5 EMA selected 0.237729 0.661055
10 TOC1 EMA selected 0.157152 1.000000
11 LUX EMA selected 0.309274 0.508131
12 ELF4 EMA selected 0.309274 0.508131
13 GI EMA selected 0.453790 0.346310
In [183]:
plt.style.use('seaborn-white')
figure(figsize=(10,6))
sns.set_palette("Paired")
pbm_vs_U2017 =  read_excel('data/PBM_vs_U20172.xlsx')
sns.barplot(x='Gene promoter', y='ratio',data=pbm_vs_U2017, hue='Model')
legend(loc='upper left', fontsize=16)
ylabel('$k_d/k_d^{TOC1}$', fontsize=22)
xlabel('Gene promoter',fontsize=22)
title('B', fontsize=22)
xticks(fontsize=16)
yticks(fontsize=16)
savefig('Kd_comparision.pdf', format='pdf', dpi=300)
In [280]:
U2017 = te.loada('Models/U2017.txt')
#
U2017.reset() #= te.loada('Models/U2017.txt')
for i in U2017p.keys():
    U2017.setValue(i,U2017p.getByKey(i))

U2017.setValue('ge3',1)
range(0,24*10,1)
res = U2017.simulate(0,24*15,24*15)
#plot(res['time'],res['[cLUX]']/max(res['[cLUX]']))
U2017 = te.loada('Models/U2017.txt')

U2017p = Utility.load('Models/params_4.bp')
U2017sf = Utility.load('Models/sfU2017.bp')
for i in U2017p.keys():
    U2017.setValue(i,U2017p.getByKey(i))
for i in U2017sf.keys():
    U2017.setValue(i,U2017sf.getByKey(i))
for i in p.keys():
    U2017.setValue(i,p[i])

res = U2017.simulate(0,24*15,3000)


figure(figsize=(15,5))

subplot(121)
plot(res['time']-24*11,res['[cG_m]'],'k--', label='cG_m')
print 'G', max(res['[cG_m]'][1500:])
#yscale('log')
xticks(range(0,24*2,12), fontsize=16)
yticks(fontsize=16)
xlim(-24,24*2)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
legend(loc='upper right', fontsize=16)
ylabel('#Transcripts/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)

ylim(0,400)
subplot(122)
plot(res['time']-24*11,res['[cT_m]'],'b--',label='cT_m')
print 'T', max(res['[cT_m]'][1500:])



U2017.setValue('g15',0.57)
U2017.setValue('g5',0.4)
res = U2017.simulate(0,24*15,3000)

subplot(121)
plot(res['time']-24*11,res['[cG_m]'],'k', label='cG_m')
print 'GI', max(res['[cG_m]'][1500:])
legend(loc='upper right', fontsize=16)
subplot(122)
plot(res['time']-24*11,res['[cT_m]'],'b',label='cT_m')
legend(loc='upper right', fontsize=16)
print 'TOC', max(res['[cT_m]'][1500:])


#yscale('log')
xticks(range(0,24*2,12), fontsize=16)
yticks(fontsize=16)
ylim(0,400)
xlim(-24,24*2)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.3)
axvspan(12,24, color='gray',alpha=0.3)
ylabel('#Transcripts/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)

#savefig('GI_and_TOC1_params_updated.png', format='png',dpi=300)

G 237.70860756406006
T 71.54658805989442
GI 292.44079777486814
TOC 132.0887656011767
Out[280]:
Text(0.5,0,'Time in constant light (h)')
In [ ]:
print 'g8\t', U2017.getValue('g8')*U2017.getValue('ge3')/(Avogadro*A)*1e9
print 'g4\t', U2017.getValue('g4')*U2017.getValue('ge3')/(Avogadro*A)*1e9
print 'g2\t', U2017.getValue('g2')*U2017.getValue('ge3')/(Avogadro*A)*1e9
print 'g14\t', U2017.getValue('g14')*U2017.getValue('ge3')/(Avogadro*A)*1e9
In [ ]:
print 'g1\t', U2017.getValue('g1')*U2017.getValue('gp7')/(Avogadro*A)*1e9
print 'g7\t', U2017.getValue('g7')*U2017.getValue('gp7')/(Avogadro*A)*1e9
print 'g13\t', U2017.getValue('g13')*U2017.getValue('gp7')/(Avogadro*A)*1e9
In [ ]:
print 'cP9m\tg8\t',U2017.getValue('ge3')*U2017.getValue('g8')/(Avogadro*A)
print 'cTm\tg4\t',U2017.getValue('ge3')*U2017.getValue('g4')/(Avogadro*A)
print 'cE4m\tg2\t',U2017.getValue('ge3')*U2017.getValue('g2')/(Avogadro*A)
print 'cGm\tg14\t',U2017.getValue('ge3')*U2017.getValue('g14')/(Avogadro*A)

Evening complex

In [ ]:
mean(array([0.195076690644,
0.360900885544,
0.158083538406,
0.104501865224,
0.20565910948,
0.223254495126,
0.301758157905]))

PBM CCA1 concentration

In [210]:
figure(figsize=(7,5))
x = linspace(0,20,100)
used = 1e-6/(42.5+66.975)*1e9
plot(x,x**2/(x**2+2.5**2), label = 'CCR2')
plot(x,x**2/(x**2+5**2), label='wt PRR9')
plot(x,x**2/(x**2+11**2), label='PRR9 CBS')
axvline(used, color='red', label='PBM concentration')
ylabel('P(CCA1 bound)', fontsize=22)
xlabel('MPB-CCA1 nM', fontsize=22)
yticks(fontsize=22)
xticks(fontsize=22)
tight_layout()
legend(loc='lower right', fontsize=16)
savefig('CCA1_PBM_affinity.png', format='png', dpi=300)
In [ ]:
figure(figsize=(15,5))
subplot(121)
AT = linspace(1e-9,1e-6,1000000)
BT = 1100/(Avogadro*1.12e-13)
k=0.35*1e-9
A = 0.5*(AT-BT-k+sqrt((AT-BT-k)**2+4*AT*k))
plot(AT,A,'blue')
plot(AT,AT,'red')
yscale('log')
xscale('log')
#ylim(1e-2,1e5)
#xlim(1e0,1e5)
xticks(fontsize=22)
yticks(fontsize=22)
ylabel('Active cL nM', fontsize=22)
xlabel('nM', fontsize=22)
subplot(122)
prr9_kd = 0.355309957403*1e-9
plot(AT,A/(prr9_kd+A), color='blue', label='Titration')
plot(AT,AT**2/(prr9_kd**2+AT**2), color='red', label='Graded')
xscale('log')
xlabel('nM')
#xlim(1,1e3)
ylim(0,1.2)
legend(loc='lower right', fontsize=16)
xticks(fontsize=22)
yticks(fontsize=22)
ylabel('cL binding probability', fontsize=22)
xlabel('nM', fontsize=22)
In [ ]:
cL = res['[cL]']/(Avogadro*1.12e-13)
plot(res['time'],cL)
In [ ]:
(U2017.getValue('gl')*U2017.getValue('g9'))/(Avogadro*A)
In [ ]:
A = 0.1168e-12
#plot(res['time'],res['[cL]']**2/(k**2+res['[cL]']**2),'o')
cL_con = 0.5*(cL-BT-k+sqrt((cL-BT-k)**2+4*cL*k))
plot(res['time'],cL_con/((U2017.getValue('gl')*U2017.getValue('g9')/(Avogadro*A)/1100+cL_con)),'-')
#yscale('log')
In [ ]:
U2017.getValue('gl')*U2017.getValue('g9')/(Avogadro*A)
In [ ]:
figure(figsize=(10,6))
xticks(range(0,24*2,12),fontsize=16)
plot(res['time']-24*11,res['[cEC]'], 'k')

yticks(fontsize=16)
xlim(-24,24*2)
#ylim(0,3e4)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.1)
axvspan(12,24, color='gray',alpha=0.1)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
title('cEC', fontsize=22)
savefig('EC_dynamics_rescaled.png', format='png',dpi=300)
In [ ]:
figure(figsize=(10,6))
xticks(range(0,24*2,12),fontsize=16)
plot(res['time']-24*11,res['[cT_m]']*1e2, 'k')
plot(res['time']-24*11,res['[cEC]'], 'k')



yticks(fontsize=16)
xlim(-24,24*2)
#ylim(0,3e4)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.1)
axvspan(12,24, color='gray',alpha=0.1)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
title('cEC', fontsize=22)
In [ ]:
figure(figsize=(10,6))
xticks(range(0,24*2,12),fontsize=16)
plot(res['time']-24*11,res['[cT_m]']*1e2, 'k')
plot(res['time']-24*11,res['[cEC]'], 'k')



yticks(fontsize=16)
xlim(-24,24*2)
#ylim(0,3e4)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.1)
axvspan(12,24, color='gray',alpha=0.1)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
title('cEC', fontsize=22)
In [ ]:
figure(figsize=(10,6))
xticks(range(0,24*2,12),fontsize=16)
U2017.reset()
res = U2017.simulate(0,24*15,3000)
plot(res['time']-24*11,res['[cT_m]']*1e2, 'k')
plot(res['time']-24*11,res['[cEC]'], 'k')



yticks(fontsize=16)
xlim(-24,24*2)
#ylim(0,3e4)
axvspan(-12,0, color='gray',alpha=0.7)
axvspan(24+12,24*2, color='gray',alpha=0.1)
axvspan(12,24, color='gray',alpha=0.1)
ylabel('#monomers/cell', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
title('cEC', fontsize=22)
In [ ]:
U2017 = te.loada('Models/U2017.txt')
#

U2017 = te.loada('Models/U2017.txt')

U2017p = Utility.load('Models/params_4.bp')
U2017sf = Utility.load('Models/sfU2017.bp')
for i in U2017p.keys():
    U2017.setValue(i,U2017p.getByKey(i))
for i in U2017sf.keys():
    U2017.setValue(i,U2017sf.getByKey(i))
for i in p.keys():
    U2017.setValue(i,p[i])

res = U2017.simulate(0,24*15,3000)
In [ ]:
cL_temp = linspace(1,1e3)
hill = U2017.p3*U2017.gl*cL_temp**2/(U2017.gl*U2017.g3**2+cL_temp**2)+cL_temp*U2017.m3
linear = [res['[cL_m]'][argmax(res['[cL]'])]*(U2017.p1+U2017.p2)/U2017.gml]*len(cL_temp)
In [ ]:
plot(cL_temp,hill)
plot(cL_temp,linear)
xscale('log')
In [ ]:
k = []
for i in pro_data['cL'].keys():
    k.append(pro_data['cL'][i][0])
k=array(k)
In [ ]:
#in aminoacids 
amin_sec = 3.
rib_tranc = 10.
CCA1_length = 609.
CCA1_trans_rate = 1/(CCA1_length/amin_sec/10)*3600
print CCA1_trans_rate
amin_sec = 3.
rib_tranc = 10.
LHY_length = 646.
LHY_trans_rate = 1/(LHY_length/amin_sec/10)*3600
print LHY_trans_rate*exp(lhyinter(sorted(pro_data['cL'].keys())))
In [1208]:
figure(figsize=(7,5))
p = []
for i in pro_data['cL'].keys():
    p.append(pro_data['cL'][i][0])

    
synth =LHY_trans_rate*exp(lhyinter(sorted(pro_data['cL'].keys())))

deg = (1.76)*array(p)

plot(array(sorted(pro_data['cL'].keys()))-24*4-24,(synth/deg), 'b',label='SM')
xlim(0,24*3)

synth = res['[cL_m]']*(U2017.p1+U2017.p2)*U2017.gl/U2017.gml
deg = U2017.m2*res['[cL]']+U2017.p3*U2017.gl*res['[cL]']**2/((U2017.gl+U2017.g3)**2+res['[cL]']**2)

plot(res['time']-24*10-24,(synth/deg), 'g',label='U2017')
xlim(0,24*4)
xticks(range(0,24*3,12),fontsize=16)
yticks(fontsize=16)
xlim(-24,24*2)
yscale('log')
axvspan(12-24,24-24,0, color='gray',alpha=0.7)
axvspan(12+24-24,24*2-24, color='gray',alpha=0.1)
axvspan(12+24*2-24,24*3-24, color='gray',alpha=0.1)
ylabel('log(synt/deg)', fontsize=22)
xlabel('Time in constant light (h)', fontsize=22)
#title('A', fontsize=22)
legend(loc='upper right', fontsize=16)
savefig('cL_synth_deg.png', format='png',dpi=300)
In [ ]:
print log10(max(res['[cL]'][1000:1300]))
print log10(min(res['[cL]'][1000:1300]))
In [ ]:
synth = res['[cL_m]']*(U2017.p1+U2017.p2)*U2017.gl/U2017.gml
deg = U2017.m2*res['[cL]']+U2017.p3*U2017.gl*res['[cL]']**2/((U2017.gl+U2017.g3)**2+res['[cL]']**2)
In [ ]:
plot(res['time']-24*10,(synth-deg))
xlim(0,24*4)
In [ ]:
In [ ]: