Curve Fitting for Reversible Two Object Binding

data_to_fit.txt

0	0.145	0.33495	0.77285	1.7835	3.6685	9.541	22.04	50.895	117.595	271.15	481.4	1450
339.3333333333	492	1418.6666666667	1341	3130	5509.3333333333	11389.3333333333	19671.6666666667	28406	35452	36769.6666666667	41462.8422425033	44774.3333333333
20.4042479237	35.9304884464	968.3849096993	125.9325216138	688.8454108144	573.4407845047	508.4302639825	691.4913834122	1150.2830086548	1451.6841943067	1052.0058618341	1456.6797863635	437.9866816849
0.1
10.0
440000.0

fitting_script.py

#!/usr/bin/env python3

#Paul Barton
#This script is in a very raw state, it could at least do with some
#improvements to inputting and to plotting.

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as pyplot


def theoretical_function(titrated_total, constant_total, kd, scale):
    a = 1.0
    b = -1.0 * (kd + titrated_total + constant_total)
    c = titrated_total * constant_total

    determinant = np.power(b, 2) - (4 * a * c)
    raw = ((-1.0 * b) - np.sqrt(determinant)) / (2.0 * a)
    return scale * raw


if __name__  == '__main__':
    with open('data_to_fit.txt', 'r') as inpt:
        titrated_values = np.array(inpt.readline().split(), dtype=float)
        measurements = np.array(inpt.readline().split(), dtype=float)
        stddevs = np.array(inpt.readline().split(), dtype=float)
        constant_value = float(inpt.readline())
        kd_guess = float(inpt.readline())
        scale_guess = float(inpt.readline())

    def fitting_theoretical(titrated_total, kd, scale):
        return theoretical_function(titrated_total,
                                    constant_value,
                                    kd,
                                    scale)

    popt, pcov = curve_fit(fitting_theoretical,  # f
                           titrated_values,  # xdata
                           measurements,  # ydata
                           p0=[kd_guess, scale_guess],
                           sigma=stddevs
                           )

    print(popt)
    
    theoretical_titration = np.linspace(titrated_values[0], titrated_values[-1], len(titrated_values) * 10)

    pyplot.title('Fitting Data to General Equation for Binding')
    pyplot.plot(titrated_values, measurements, 'rd', label='Measured Data')
    pyplot.plot(theoretical_titration, fitting_theoretical(theoretical_titration, *popt), 'b-', label='Fit Curve')
    pyplot.xlabel('Titrated Agent Concentration (uM)')
    pyplot.ylabel('Relative Fluorescence Units (RFU)')
    pyplot.legend(loc='best')
    pyplot.savefig('plot.png')