In this notebook we will see how well we can reproduce Kd of a non-fluorescent ligand from simulated experimental data with a maximum likelihood function.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
import seaborn as sns
%pylab inline
In [2]:
#Competitive binding function
#This function and its assumptions are defined in greater detail in this notebook:
## modelling-CompetitiveBinding-ThreeComponentBinding.ipynb
def three_component_competitive_binding(Ptot, Ltot, Kd_L, Atot, Kd_A):
"""
Parameters
----------
Ptot : float
Total protein concentration
Ltot : float
Total tracer(fluorescent) ligand concentration
Kd_L : float
Dissociation constant of the fluorescent ligand
Atot : float
Total competitive ligand concentration
Kd_A : float
Dissociation constant of the competitive ligand
Returns
-------
P : float
Free protein concentration
L : float
Free ligand concentration
A : float
Free ligand concentration
PL : float
Complex concentration
Kd_L_app : float
Apparent dissociation constant of L in the presence of A
Usage
-----
[P, L, A, PL, Kd_L_app] = three_component_competitive_binding(Ptot, Ltot, Kd_L, Atot, Kd_A)
"""
Kd_L_app = Kd_L*(1+Atot/Kd_A)
PL = 0.5 * ((Ptot + Ltot + Kd_L_app) - np.sqrt((Ptot + Ltot + Kd_L_app)**2 - 4*Ptot*Ltot)) # complex concentration (uM)
P = Ptot - PL; # free protein concentration in sample cell after n injections (uM)
L = Ltot - PL; # free tracer ligand concentration in sample cell after n injections (uM)
A = Atot - PL; # free competitive ligand concentration in sample cell after n injections (uM)
return [P, L, A, PL, Kd_L_app]
In [3]:
#Let's define our parameters
Kd = 3800e-9 # M
Kd_Competitor = 3000e-9 # M
Ptot = 1e-9 * np.ones([12],np.float64) # M
Ltot = 20.0e-6 / np.array([10**(float(i)/2.0) for i in range(12)]) # M
L_Competitor = 10e-6 # M
In [4]:
[P, L, A, PL, Kd_L_app] = three_component_competitive_binding(Ptot, Ltot, Kd, L_Competitor, Kd_Competitor)
#usin
[P_base, L_base, A_base, PL_base, Kd_L_app_base] = three_component_competitive_binding(Ptot, Ltot, Kd, 0, Kd_Competitor)
In [5]:
# y will be complex concentration
# x will be total ligand concentration
plt.title('Competition assay')
plt.semilogx(Ltot,PL_base,'green', marker='o', label = 'Fluorescent Ligand')
plt.semilogx(Ltot,PL,'cyan', marker='o', label = 'Fluorescent Ligand + Competitor')
plt.xlabel('$[L]_{tot}$')
plt.ylabel('$[PL]$')
plt.legend(loc=0);
In [6]:
#What if we change our Kd's a little
#Kd_Gef_Abl
Kd = 480e-9 # M
#Kd_Ima_Abl
Kd_Competitor = 21.0e-9 # M
In [7]:
[P, L, A, PL, Kd_L_app] = three_component_competitive_binding(Ptot, Ltot, Kd, L_Competitor, Kd_Competitor)
#using _base as a subscript to define when we have no competitive ligand
[P_base, L_base, A_base, PL_base, Kd_L_app_base] = three_component_competitive_binding(Ptot, Ltot, Kd, 0, Kd_Competitor)
In [8]:
# y will be complex concentration
# x will be total ligand concentration
plt.title('Competition assay')
plt.semilogx(Ltot,PL_base,'green', marker='o', label = 'Fluorescent Ligand')
plt.semilogx(Ltot,PL,'cyan', marker='o', label = 'Fluorescent Ligand + Competitor')
plt.xlabel('$[L]_{tot}$')
plt.ylabel('$[PL]$')
plt.legend(loc=0);
In [9]:
# Making max 400 relative fluorescence units, and scaling all of PL to that
npoints = len(Ltot)
sigma = 10.0 # size of noise
F_i = (400/1e-9)*PL + sigma * np.random.randn(npoints)
F_i_base = (400/1e-9)*PL_base + sigma * np.random.randn(npoints)
In [10]:
# y will be complex concentration
# x will be total ligand concentration
plt.title('Competition assay')
plt.semilogx(Ltot,F_i_base,'green', marker='o', label = 'Fluorescent Ligand')
plt.semilogx(Ltot,F_i,'cyan', marker='o', label = 'Fluorescent Ligand + Competitor')
plt.xlabel('$[L]_{tot}$')
plt.ylabel('$Fluorescence$')
plt.legend(loc=0);
In [11]:
#And makeup an F_L
F_L = 0.3
In [12]:
# This function fits Kd_L when L_Competitor is 0
def find_Kd_from_fluorescence_base(params):
[F_background, F_PL, Kd_L] = params
N = len(Ltot)
Fmodel_i = np.zeros([N])
for i in range(N):
[P, L, A, PL, Kd_L_app] = three_component_competitive_binding(Ptot[0], Ltot[i], Kd_L, 0, Kd_Competitor)
Fmodel_i[i] = (F_PL*PL + F_L*L) + F_background
return Fmodel_i
In [ ]:
In [13]:
initial_guess = [1,400/1e-9,3800e-9]
prediction = find_Kd_from_fluorescence_base(initial_guess)
plt.semilogx(Ltot,prediction,color='k')
plt.semilogx(Ltot,F_i_base, 'o')
plt.xlabel('$[L]_{tot}$ / M')
plt.ylabel('$Fluorescence$')
plt.legend();
In [14]:
def sumofsquares(params):
prediction = find_Kd_from_fluorescence_base(params)
return np.sum((prediction - F_i_base)**2)
In [15]:
initial_guess = [0,3E11,2000E-9]
fit = optimize.minimize(sumofsquares,initial_guess,method='Nelder-Mead')
print "The predicted parameters are", fit.x
In [16]:
fit_prediction = find_Kd_from_fluorescence_base(fit.x)
In [17]:
plt.semilogx(Ltot,fit_prediction,color='k')
plt.semilogx(Ltot,F_i_base, 'o')
plt.xlabel('$[L]_{tot}$ / M')
plt.ylabel('$Fluorescence$')
plt.legend();
In [18]:
Kd_L_MLE = fit.x[2]
In [19]:
def Kd_format(Kd):
if (Kd < 1e-12):
Kd_summary = "Kd = %.1f nM " % (Kd/1e-15)
elif (Kd < 1e-9):
Kd_summary = "Kd = %.1f pM " % (Kd/1e-12)
elif (Kd < 1e-6):
Kd_summary = "Kd = %.1f nM " % (Kd/1e-9)
elif (Kd < 1e-3):
Kd_summary = "Kd = %.1f uM " % (Kd/1e-6)
elif (Kd < 1):
Kd_summary = "Kd = %.1f mM " % (Kd/1e-3)
else:
Kd_summary = "Kd = %.3e M " % (Kd)
return Kd_summary
In [20]:
Kd_format(Kd_L_MLE)
Out[20]:
In [21]:
delG_summary = "delG = %s kT" %np.log(Kd_L_MLE)
delG_summary
Out[21]:
In [22]:
# This function fits Kd_A when Kd_L already has an estimate
def find_Kd_from_fluorescence_competitor(params):
[F_background, F_PL, Kd_Competitor] = params
N = len(Ltot)
Fmodel_i = np.zeros([N])
for i in range(N):
[P, L, A, PL, Kd_L_app] = three_component_competitive_binding(Ptot[0], Ltot[i], Kd_L_MLE, L_Competitor, Kd_Competitor)
Fmodel_i[i] = (F_PL*PL + F_L*L) + F_background
return Fmodel_i
In [23]:
initial_guess = [0,400/1e-9,3800e-9]
prediction = find_Kd_from_fluorescence_competitor(initial_guess)
plt.semilogx(Ltot,prediction,color='k')
plt.semilogx(Ltot,F_i, 'o')
plt.xlabel('$[L]_{tot}$ / M')
plt.ylabel('$Fluorescence$')
plt.legend();
In [24]:
def sumofsquares(params):
prediction = find_Kd_from_fluorescence_competitor(params)
return np.sum((prediction - F_i)**2)
In [25]:
initial_guess = [0,3E11,2000E-9]
fit_comp = optimize.minimize(sumofsquares,initial_guess,method='Nelder-Mead')
print "The predicted parameters are", fit_comp.x
In [26]:
fit_prediction_comp = find_Kd_from_fluorescence_competitor(fit_comp.x)
In [27]:
plt.semilogx(Ltot,fit_prediction_comp,color='k')
plt.semilogx(Ltot,F_i, 'o')
plt.xlabel('$[L]_{tot}$ / M')
plt.ylabel('$Fluorescence$')
plt.legend();
In [28]:
Kd_A_MLE = fit_comp.x[2]
In [29]:
Kd_format(Kd_A_MLE)
Out[29]:
In [30]:
delG_summary = "delG = %s kT" %np.log(Kd_A_MLE)
delG_summary
Out[30]:
In [31]:
plt.semilogx(Ltot,fit_prediction_comp,color='k')
plt.semilogx(Ltot,F_i, 'o')
plt.axvline(Kd_A_MLE,color='k',label='%s (MLE)'%Kd_format(Kd_A_MLE))
plt.axvline(Kd_Competitor,color='b',label='%s'%Kd_format(Kd_Competitor))
plt.semilogx(Ltot,fit_prediction,color='k')
plt.semilogx(Ltot,F_i_base, 'o')
plt.axvline(Kd_L_MLE,color='k',label='%s (MLE)'%Kd_format(Kd_L_MLE))
plt.axvline(Kd,color='g',label='%s'%Kd_format(Kd))
plt.xlabel('$[L]_{tot}$ / M')
plt.ylabel('$Fluorescence$')
plt.legend(loc=0);
In [32]:
#Awesome
In [ ]: