In [1]:
%pylab inline

import numpy as np
import scipy
import math


Populating the interactive namespace from numpy and matplotlib

In [2]:
def calc_abc(kd, p_total, l_total, alpha):
    a = 2.0 * kd / alpha + 2.0 * l_total - p_total
    b = (math.pow(kd, 2.0) + 2 * kd * l_total - 2 * kd * p_total) / alpha
    c = -1 * (math.pow(kd, 2) * p_total) / alpha
    return a, b, c

In [3]:
def calc_qr(a, b, c):
    q = (3 * b - math.pow(a, 2.0)) / 9
    r = (9 * a * b - 27 * c - 2 * math.pow(a, 3.0)) / 54
    return q, r

In [4]:
def cartesian_cubic(a, q, r):
    first = -1 * a / 3.0
    second = math.pow(r + math.pow(math.pow(q, 3.0) + math.pow(r, 2.0), 0.5), 1.0 / 3.0)
    third = math.pow(r - math.pow(math.pow(q, 3.0) + math.pow(r, 2.0), 0.5), 1.0 / 3.0)
    return first + second + third

def polar_cubic(a, q, r):
    theta =  math.acos(r / math.pow(-1 * math.pow(q, 3), 1.0 / 2.0))
    return math.cos(theta / 3.0) * math.pow(-1 * q, 1.0 / 2.0) * 2 - (a / 3.0)

In [5]:
def model(kd, p_total, l_total, alpha):
    a, b, c = calc_abc(kd, p_total, l_total, alpha)
    q, r = calc_qr(a, b, c)
    #Use of cartesian or cubic depends on Q^3 + R^2
    if math.pow(q, 3.0) + math.pow(r, 2.0) > 0:
         return cartesian_cubic(a,q,r)
    else:
        print(polar_cubic(a, q, r))
        return polar_cubic(a, q, r)

In [6]:
#Normalize all values to uM
kd = 1
p_total = .1
alpha = 1

#Create the plot
plot = pylab.figure().add_subplot(111)

lig_range = [0.00001 * math.pow(10, lig_val) for lig_val in np.linspace(1, 10, 91)]
model_vals = [model(kd, p_total, l_total, alpha) / p_total for l_total in lig_range]

plot.plot(lig_range, model_vals)
plot.set_ylabel('[P]')
plot.set_xscale('log')
plot.set_xlabel(r'Ligand (uM)')
plot.grid()


0.09998182118663046
0.09997711520918129
0.09997119130722754
0.09996373446396634
0.09995434829242944
0.09994253408100284
0.09992766448003654
0.09990895048088444
0.09988540001214075
0.09985576607980173
0.09981848189762121
0.09977157988152197
0.09971259071223315
0.09963841790694927
0.09954518250412447
0.0994280316053412
0.09928090373368315
0.09909624343804124
0.09886465760107177
0.09857450698834214
0.09821142947811456
0.09775779729340806
0.09719212108600284
0.0964884311580596
0.09561569327544162
0.09453735647029071
0.09321118522815441
0.09158959790541576
0.08962080979873699
0.08725114305336645
0.08442887702247603
0.08110990857711675
0.07726519341610671
0.07288938669152634
0.06800932439958196
0.06269019700191225
0.05703689131967171
0.05118852097551296
0.0453058488212974
0.039553674071090006
0.034082207965583144
0.02901188945575517
0.024424729746281626
0.02036291790958078
0.016833299039193506
0.013815263923823817
0.011269620792015989
0.009146711193775747
0.007392873694965552
0.00595503245647766
0.00478359412570839
0.0038340114952575277
0.0030673970397909756
0.0024505181080449745
0.0019554285842708907
0.0015589166323941583
0.0012418861987804064
0.000988743790109936
0.0007868300035269726
0.000625914332964328
0.000497758846691454
0.00039574890126914397
0.0003145851555785839
0.0002500295573781841
0.00019869754484602709
0.00015788916641668038
0.00012545250280027176
9.967361017970688e-05
7.918820045915709e-05
6.291058537044592e-05
4.997747510060435e-05
3.970223326632549e-05
3.153888519591419e-05
2.5053930357898935e-05
1.9902170379282325e-05
1.5809526303200983e-05
1.2557242371258326e-05
9.976409728551516e-06
7.925065801828168e-06
6.291312274697702e-06
5.014466296415776e-06
3.975330400862731e-06
3.1461240723729134e-06
2.510108970454894e-06
1.947966666193679e-06
1.6051089914981276e-06
1.1585289030335844e-06
1.0417134035378695e-06
5.966285243630409e-07
5.377587513066828e-07
-2.0960578694939613e-07

In [ ]: