initial work on the data fitting

data_fitting/fig7/f7minus_auto.csv

+,Image Name,Channel,Time,Signal,Total,Area,Bkgnd.,Type
+0,070119 shka d2a 1.5 hrs.tif,-cdG,15,35386532.4250115,44908956,7056,1349.54968261719,Signal
+1,070119 shka d2a 1.5 hrs.tif,,30,52905714.4250115,62428136,7056,1349.54968261719,Signal
+2,070119 shka d2a 1.5 hrs.tif,,60,74998865.4250115,84521288,7056,1349.54968261719,Signal
+3,070119 shka d2a 1.5 hrs.tif,,90,87397755.4250115,96920176,7056,1349.54968261719,Signal
+4,070119 shka d2a 1.5 hrs.tif,,180,105658743.42501199,115181168,7056,1349.54968261719,Signal
+5,070119 shka d2a 1.5 hrs.tif,,360,125704783.425012,135227200,7056,1349.54968261719,Signal
+6,070119 shka d2a 1.5 hrs.tif,,720,130330163.425012,139852592,7056,1349.54968261719,Signal

data_fitting/fig7/f7minus_fplc.csv

+,run,time [sec],E [M],S,S [M],remarks,cv(??),cv_fit(??),area(??),cv(Prot),cv_fit(Prot),area(Prot),cv(ADP),cv_fit(ADP),area(ADP),cv(ATP),cv_fit(ATP),area(ATP),cv(cdG),cv_fit(cdG),area(cdG),sum(nucl)
+0,2,57,1e-05,ATP,0.0002,,2,2.05,-0.3091237798875516,3.6,3.55,2.56492752158952,4.52,4.530423103103309,0.971082217952269,5.31,5.29608587648113,66.9145704657649,7.33,7.38,-0.2294294543330736,67.8856526837172
+0,1,90,1e-05,ATP,0.0002,,2,2.0,-0.2178417057408677,3.55,3.5,2.6989074665730106,4.57,4.52,1.239649690915134,5.26,5.2279550365085,64.2037447678706,7.33,7.38,-0.2399302413163968,65.4433944587858
+1,3,318,1e-05,ATP,0.0002,,2,2.05,-0.2943369482397821,3.6,3.56699896802388,2.723785446126083,4.52,4.510186703449049,1.602731245702074,5.31,5.291972420416941,66.0149234087844,7.33,7.38,-0.2732121351178214,67.6176546544865
+1,2,350,1e-05,ATP,0.0002,,2,2.05,-0.330586136015984,3.55970067129851,3.55970219615754,2.7326041005108146,4.50321962459605,4.503220241196059,1.656565870508687,5.28679968792753,5.2867997354002805,64.4060393416747,7.28,7.255384946880669,-0.2569803203715678,66.0626052121834
+2,4,579,1e-05,ATP,0.0002,,2,2.05,-0.328000781969354,3.6,3.5686517687111903,2.733741474553166,4.52,4.50780518970171,1.6487861355283668,5.31,5.2910964227249195,65.3126535012086,7.33,7.28,-0.2851824186194056,66.961439636737
+2,3,610,1e-05,ATP,0.0002,,2,2.0,-0.1926287255753005,3.55,3.5,2.7814357383574158,4.47,4.43695013323787,1.814050637136272,5.26,5.2192633803520305,65.3016722140969,7.33,7.38,-0.32314799289930396,67.1157228512332
+3,5,901,1e-05,ATP,0.0002,,2,2.0,-0.327126290854358,3.6,3.55,2.648228994431745,4.52,4.47,1.7535223922980332,5.31,5.26,65.643533345465,7.33,7.33,-0.0985258318967615,67.397055737763
+3,4,930,1e-05,ATP,0.0002,,2,2.05,-0.2435775013799569,3.56335071375361,3.5633476510947504,2.727069171746906,4.501192380008179,4.50118737069838,1.774537069596054,5.28265861518206,5.282658650329941,66.3081402195577,7.28000214903484,7.258658272487481,-0.25507563070821,68.0826772891537
+4,6,1343,1e-05,ATP,0.0002,,2,2.05,-0.31641803160173404,3.6,3.5726595269724704,2.7271012187832433,4.52,4.511346098031191,1.7540172561046108,5.31,5.29226782023431,65.2715025489618,7.33,7.38,-0.2509514469830947,67.0255198050664
+4,5,1370,1e-05,ATP,0.0002,,2,2.05,-0.2959208327889744,3.55,3.50470261760739,2.736272770509535,4.47,4.44353203377238,1.874997716200741,5.26,5.22642306013939,67.242617673452,7.33,7.38,-0.2914350294277524,69.1176153896528

data_fitting/fig7/f7plus_auto.csv

+,Image Name,Channel,Time,Signal,Total,Area,Bkgnd.,Type
+0,070119 shka d2a 1.5 hrs.tif,+cdG,15,175063445.42501202,184585872,7056,1349.54968261719,Signal
+1,070119 shka d2a 1.5 hrs.tif,,30,181952993.42501202,191475424,7056,1349.54968261719,Signal
+2,070119 shka d2a 1.5 hrs.tif,,60,185346201.42501202,194868624,7056,1349.54968261719,Signal
+3,070119 shka d2a 1.5 hrs.tif,,90,180689021.42501202,190211440,7056,1349.54968261719,Signal
+4,070119 shka d2a 1.5 hrs.tif,,180,182724421.42501202,192246848,7056,1349.54968261719,Signal
+5,070119 shka d2a 1.5 hrs.tif,,360,192085249.42501202,201607680,7056,1349.54968261719,Signal
+6,070119 shka d2a 1.5 hrs.tif,,720,195943527.42501202,205465952,7056,1349.54968261719,Signal

data_fitting/fig7/f7plus_fplc.csv

+,run,time [sec],E [M],S,S [M],remarks,cv(??),cv_fit(??),area(??),cv(Prot),cv_fit(Prot),area(Prot),cv(ADP),cv_fit(ADP),area(ADP),cv(ATP),cv_fit(ATP),area(ATP),cv(cdG),cv_fit(cdG),area(cdG),sum(nucl)
+0,2,56,1e-05,ATP,0.0002,,2,2.05,-0.40112313692673296,3.6,3.60232471358526,2.771955094337488,4.52,4.50924404612424,2.2011749700078433,5.31,5.291856500579049,65.3775787777451,7.33,7.32696295147421,6.92977461244807,67.578753747753
+0,1,90,1e-05,ATP,0.0002,,2,2.05,-0.335549784033709,3.6,3.5992858379486203,2.728961040132289,4.52,4.51176375436566,2.222913823341912,5.31,5.295223810650549,64.35885425420629,7.33,7.33326815337582,7.4429749301643895,66.5817680775482
+1,3,318,1e-05,ATP,0.0002,,2,2.05,-0.39969911264739,3.6,3.606193084799,2.7541620850407673,4.52,4.51120128585202,2.2690117673785837,5.31,5.29253329041674,64.8778124645464,7.33,7.329689962327611,7.20969383578719,67.146824231925
+1,2,350,1e-05,ATP,0.0002,,2,2.05,-0.42804789260648207,3.6,3.6156661783619906,2.7624177568276966,4.52,4.52442454407732,2.323452091451121,5.31,5.31552675107883,63.809764920275796,7.33,7.3430373654480805,7.680659483689059,66.1332170117269
+2,4,579,1e-05,ATP,0.0002,,2,2.05,-0.39007028537106603,3.6,3.5979159652895403,2.821655528222326,4.52,4.50299379466243,2.315543652051771,5.31,5.28613708974721,64.65328354440621,7.33,7.32093722426105,7.299559160902151,66.968827196458
+2,3,610,1e-05,ATP,0.0002,,2,2.05,-0.39710935060842895,3.6,3.6097891028273397,2.6886333423169577,4.52,4.52474384999891,2.2546417781887405,5.31,5.31613655704751,65.0603868484838,7.33,7.34777861825318,7.78250720389114,67.3150286266725
+3,5,901,1e-05,ATP,0.0002,,2,2.05,-0.42168659906366707,3.6,3.60905217367437,2.832577279673719,4.52,4.5147651880821105,2.3782888588942908,5.31,5.29860124543565,64.3924118307112,7.33,7.33234529486231,7.36209317957299,66.77070068960539
+3,4,930,1e-05,ATP,0.0002,,2,2.05,-0.3904818405954771,3.6,3.61172909674684,2.725426429919315,4.52,4.52304937927825,2.3287182270715068,5.31,5.31590747000612,65.5254659026753,7.33,7.35218334266989,7.73936108734887,67.8541841297468
+4,6,1343,1e-05,ATP,0.0002,,2,2.05,-0.430130175854087,3.6,3.61228682283454,2.8290684438565665,4.52,4.51730786108753,2.334664930475081,5.31,5.300506100367509,63.8247068575288,7.33,7.330908122191531,7.48833432417588,66.1593717880038
+4,5,1370,1e-05,ATP,0.0002,,2,2.05,-0.4044442184441911,3.6,3.6115844875937,2.7217181067553944,4.52,4.52196861866195,2.3137062650146545,5.31,5.31084964923812,66.2809846359552,7.33,7.34745333614643,7.8427135303497195,68.5946909009699

data_fitting/kinetic.py

+../kinetic.py
\ No newline at end of file

data_fitting/refine_fig4.py

+#!/struct/soft/osx/anaconda3/bin/python
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+from shka_model_monomer import concentrations_model, equilibrate_model
+
+#%%
+df_fplc  = pd.read_csv('fig7/f7minus_fplc.csv')
+df_fplc2 = pd.read_csv('fig7/f7plus_fplc.csv')
+
+df_auto  = pd.read_csv('fig7/f7minus_auto.csv')
+df_auto2 = pd.read_csv('fig7/f7plus_auto.csv')
+
+#%% normalize auto
+df_auto.Signal /= 1e8
+df_auto2.Signal /= 1e8
+
+#%% 7a
+plt.figure(figsize=(10,5))
+plt.subplot(1,2,1)
+plt.plot(df_auto.Time, df_auto.Signal, 'bx')
+plt.plot(df_auto2.Time, df_auto2.Signal, 'gx')
+#plt.xlim((0,1000))
+plt.ylim((0,3))
+
+plt.subplot(1,2,2)
+# 7c
+plt.plot(df_fplc['time [sec]'], df_fplc['area(ADP)'], 'bx')
+plt.plot(df_fplc2['time [sec]'], df_fplc2['area(ADP)'], 'rx')
+#plt.xlim((0,1000))
+plt.ylim((0,3))
+
+#%% User input:
+K1 = 0.33
+K2 = 2.0
+K01 = 50.0
+
+# initial conditions
+c0 = {'_HD.free': 5.0, 'ATP': 200.0}
+
+# rate constants
+k01  = 100.0       #   _HD -> DH=HD
+k01r = k01/K01     #   =HD -> DH_HD
+k02  = 0.0         #   ATP -> ADP + P
+
+k1    = 0.01   #   H.ATP -> Hp.ADP
+k1r   = k1/K1  #  Hp.ADP -> H.ATP
+k1non = 0.0    #      Hp -> H + P
+k2    = 0.0    # HpD.ATP -> HDp.ATP
+k2r   = k2/K2
+k3    = 0.0    # HDp.ADP -> HD.ADP + P
+k4    = 0.0    #      Dp -> D + P
+
+#%% equilibrate system
+
+Kd = 50 # uM
+kon = 100.0 # uM-1s-1 10^8 M-1s-1  diffusion limit
+koff = Kd*kon
+
+rates = { 'kon': kon,
+         'koff': koff,
+          'k01': k01,
+         'k01r': k01r}
+
+species, c1 = equilibrate_model(c0, rates)
+
+#%%
+rates = { 'kon': kon,
+         'koff': koff,
+          'k01': k01,
+         'k01r': k01r,
+          'k02': k02,
+           'k1': k1,
+          'k1r': k1r,
+        'k1non': k1non,
+           'k2': k2,
+          'k2r': k2r,
+           'k3': k3,
+           'k4': k4}
+
+#%% error function
+def errFunction(p):
+    global K1
+    global rates, df_auto, df_fplc, c1
+    global species, rxns_d1, ct1, ct1_adp, ct1_atp, ct1_boundp
+    global ct2, ct2_adp, ct2_atp, ct2_boundp
+    global t1, t2
+
+    k1, sc_auto, sc_fplc = p
+    print('Params: ', p, end='')
+
+    rates['k1'] = k1
+    rates['k1r'] = k1 / K1
+
+    t1 = np.insert(df_auto.Time.values, 0, 0)
+    species, rxns_d1, ct1, ct1_adp, ct1_atp, ct1_boundp = concentrations_model(c1, rates, t1)
+
+    t2 = np.insert(df_fplc['time [sec]'].values, 0, 0)
+    species, rxns_d1, ct2, ct2_adp, ct2_atp, ct2_boundp = concentrations_model(c1, rates, t2)
+
+    d1 = df_auto.Signal.values - sc_auto*ct1_boundp[1:]
+    d2 = df_fplc['area(ADP)'].values - sc_fplc*ct2_adp[1:]
+
+    d = np.append(d1, d2)
+
+    print('  chi^2: %8.3f'%np.sum(d**2))
+
+    return d
+
+#%%
+p0 = [0.3, 0.3, 0.4] # k1, sc_auto, sc_fplc
+
+d = errFunction(p0)
+
+#%% 7a
+_, sc_auto, sc_fplc = p0
+
+plt.figure(figsize=(10,5))
+plt.subplot(1,2,1)
+# exp
+plt.plot(df_auto.Time, df_auto.Signal, 'bx', label='-cdG')
+plt.plot(df_auto2.Time, df_auto2.Signal, 'gx', label='+cdG')
+# calc
+plt.plot(t1, sc_auto*ct1_boundp, '.-', label='bound P')
+plt.legend()
+
+#plt.xlim((0,1000))
+plt.ylim((0,3))
+
+# 7c
+plt.subplot(1,2,2)
+plt.plot(df_fplc['time [sec]'], df_fplc['area(ADP)'], 'bx', label='-cdG')
+plt.plot(df_fplc2['time [sec]'], df_fplc2['area(ADP)'], 'rx', label='+cdG')
+# calc
+plt.plot(t2, sc_fplc*ct2_adp, '.-', label='ADP')
+#plt.xlim((0,1000))
+plt.ylim((0,3))
+
+plt.legend()
+
+#%%
+from scipy.optimize import least_squares
+
+#%%
+
+res = least_squares(errFunction, p0, bounds=(0.0, np.inf))
+
+#%% optimized parameter values
+p1 = res.x
+
+#%%
+d = errFunction(p1)
+
+#%% 7a
+_, sc_auto, sc_fplc = p1
+
+plt.figure(figsize=(10,5))
+plt.subplot(1,2,1)
+# exp
+plt.plot(df_auto.Time, df_auto.Signal, 'bx', label='-cdG')
+plt.plot(df_auto2.Time, df_auto2.Signal, 'gx', label='+cdG')
+# calc
+plt.plot(t1, sc_auto*ct1_boundp, '.-', label='bound P')
+plt.legend()
+
+#plt.xlim((0,1000))
+plt.ylim((0,3))
+
+# 7c
+plt.subplot(1,2,2)
+plt.plot(df_fplc['time [sec]'], df_fplc['area(ADP)'], 'bx', label='-cdG')
+plt.plot(df_fplc2['time [sec]'], df_fplc2['area(ADP)'], 'rx', label='+cdG')
+# calc
+plt.plot(t2, sc_fplc*ct2_adp, '.-', label='ADP')
+#plt.xlim((0,1000))
+plt.ylim((0,3))
+
+plt.legend()

data_fitting/shka_model_monomer.py

+# -*- coding: utf-8 -*-
+
+import re
+import numpy as np
+
+import kinetic
+
+#%%
+
+def concentrations_model(c0, rates, t, hasAMPPNP=False):
+    """
+    Calculate transient concentrations
+    """
+    roots = [  '_HD',  '_HpD',  '_HDp',  '_HpDp']
+
+    roots_locked = [s.replace('_', '=') for s in roots]
+
+    ligands = ['ATP', 'ADP']
+
+    if hasAMPPNP:
+        ligands.append('AMPPNP')
+
+    # ligand equilibria
+    ligand_eq_rxns_s = []
+
+    for root in roots + roots_locked:
+        for l in ligands:
+            ligand_eq_rxns_s.append('%s.%s <=> %s.free + %s    ; kon koff'%(root, l, root, l))
+
+    # active <=> locked equilibria
+    active_locked_eq_rxns_s = []
+
+    for active, locked in zip(roots, roots_locked):
+        #print(active, ' ', locked)
+        for l in ligands + ['free']:
+            active_locked_eq_rxns_s.append('%s.%s <=> %s.%s    ; k01 k01r'%(active, l, locked, l))
+
+    #
+    eq_rxns= [kinetic.reaction(rxn_s) for rxn_s in ligand_eq_rxns_s ]
+    active_species = {k:v for k,v in kinetic.get_species(eq_rxns).items() if k.count('_')>0 }
+    eq_species = kinetic.get_species(eq_rxns)
+
+    # add reactions with k1
+    k1_rxns_s = []
+
+    k1_rxns_s.append('_HD.ATP <=> _HpD.ADP   ; k1 k1r')
+    k1_rxns_s.append('_HDp.ATP <=> _HpDp.ADP   ; k1 k1r')
+
+    # non-specific Hp dephosphorylation
+    k1non_rxns_s = []
+
+    for spc, spc_id in eq_species.items():
+        idxs = [match.start() for match in re.finditer('Hp', spc)]
+        for idx in idxs:
+            prefix = spc[:idx]
+            suffix = spc[idx+2:]
+            k1non_rxns_s.append('%s -> %sH%s + P    ; k1non'%(spc, prefix, suffix))
+
+    # add reactions with k2
+    k2_rxns_s = []
+
+    for spc, spc_id in active_species.items():
+        idx = spc.find('_HpD.')
+        if idx >= 0:
+            prefix = spc[:idx]
+            suffix = spc[idx+5:]
+            k2_rxns_s.append('%s <=> %s_HDp.%s    ; k2 k2r'%(spc, prefix, suffix))
+
+    # add reactions with k3
+    k3_rxns_s = []
+
+    # k3_rxns_s.append('_HDp.ADP -> _HD.ADP + P   ; k3adp')
+
+    for spc, spc_id in eq_species.items():
+        idxs = [match.start() for match in re.finditer('_HDp.', spc)]
+        for idx in idxs:
+            suffix = spc[idx+5:]
+            # k3_rxns_s.append('%s -> _HD.%s + P    ; k3any'%(spc, suffix))
+            k3_rxns_s.append('%s -> _HD.%s + P    ; k3'%(spc, suffix))
+
+#    for spc, spc_id in active_species.items():
+#        idx = spc.find('_HDp.ADP')
+#        if idx > 0:
+#            prefix = spc[:idx]
+#            k3_rxns_s.append('%s -> %s_HD.ADP + P ; k3'%(spc, prefix))
+#
+#        if spc.find('ADP.DpH_') >= 0:
+#            suffix = spc[8:]
+#            k3_rxns_s.append('%s -> ADP.DH_%s + P ; k3'%(spc, suffix))
+
+    # non-specific Dp dephosphorylation
+    k4_rxns_s = []
+
+    for spc, spc_id in eq_species.items():
+        idxs = [match.start() for match in re.finditer('Dp', spc)]
+        for idx in idxs:
+            prefix = spc[:idx]
+            suffix = spc[idx+2:]
+            k4_rxns_s.append('%s -> %sD%s + P    ; k4'%(spc, prefix, suffix))
+
+    # ATP hydrolysis
+    unspecific = ['ATP -> ADP + P    ; k02']
+
+    #
+    rxns = [kinetic.reaction(rxn_s) for rxn_s in ligand_eq_rxns_s + active_locked_eq_rxns_s + k1_rxns_s + k1non_rxns_s + k2_rxns_s + k3_rxns_s + k4_rxns_s + unspecific]
+
+    rxns_dict = {'ligand_eq': ligand_eq_rxns_s,
+                 'active_locked_eq': active_locked_eq_rxns_s,
+                 'k1': k1_rxns_s,
+                 'k1non': k1non_rxns_s,
+                 'k2': k2_rxns_s,
+                 'k3': k3_rxns_s,
+                 'k4': k4_rxns_s,
+                 'unspecific': unspecific}
+
+    #
+    species, ct = kinetic.concentration(rxns, c0, rates, t)
+
+    # calculate transient total ADP concentration
+    ctADP_total = np.zeros(ct.shape[0])
+    ctATP_total = np.zeros(ct.shape[0])
+    ctBoundP_total = np.zeros(ct.shape[0])
+
+    for spc, spc_id in species.items():
+        nADP = spc.count('ADP')
+        if nADP > 0:
+            ctADP_total += nADP * ct[:, spc_id]
+
+        nATP = spc.count('ATP')
+        if nATP > 0:
+            ctATP_total += nATP * ct[:, spc_id]
+
+        n_p = spc.count('p')
+        if n_p > 0:
+            ctBoundP_total += n_p * ct[:, spc_id]
+
+    return species, rxns_dict, ct, ctADP_total, ctATP_total, ctBoundP_total
+
+
+def equilibrate_model(c0, rates, dt=1.0, npoints=10, eps=1e-3, hasAMPPNP=False):
+    """
+    Equilibrate ligands binding and locked<->active states
+    """
+
+    roots = [  '_HD',  '_HpD',  '_HDp',  '_HpDp']
+
+    roots_locked = [s.replace('_', '=') for s in roots]
+
+    ligands = ['ATP', 'ADP']
+
+    if hasAMPPNP:
+        ligands.append('AMPPNP')
+
+    # ligand equilibria
+    ligand_eq_rxns_s = []
+
+    for root in roots + roots_locked:
+        for l in ligands:
+            ligand_eq_rxns_s.append('%s.%s <=> %s.free + %s    ; kon koff'%(root, l, root, l))
+
+    # active <=> locked equilibria
+    active_locked_eq_rxns_s = []
+
+    for active, locked in zip(roots, roots_locked):
+        #print(active, ' ', locked)
+        for l in ligands + ['free']:
+            active_locked_eq_rxns_s.append('%s.%s <=> %s.%s    ; k01 k01r'%(active, l, locked, l))
+
+    #
+    rxns = [kinetic.reaction(rxn_s) for rxn_s in ligand_eq_rxns_s + active_locked_eq_rxns_s]
+
+    #
+    t = np.linspace(0, dt, npoints)
+
+    species, ct = kinetic.concentration(rxns, c0, rates, t)
+
+    i = 1
+    chi2 = np.sqrt(np.sum((ct[-1,:] - ct[-2,:])**2))
+    c1 = kinetic.make_conc_dict(species, ct[-1, :])
+
+    while chi2 > eps:
+        print('iteration %d:   %f'%(i, chi2))
+        c1 = kinetic.make_conc_dict(species, ct[-1, :])
+        species, ct = kinetic.concentration(rxns, c1, rates, t)
+        chi2 = np.sqrt(np.sum((ct[-1,:] - ct[-2,:])**2))
+        i+=1
+
+    return species, c1