Some general settings:
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
import sys, time, math
import numpy as np
from dcprogs.likelihood import inv
HJCFIT depends on DCPROGS/DCPYPS module for data input and setting kinetic mechanism:
In [3]:
from dcpyps.samples import samples
from dcpyps import dataset, mechanism, dcplots
In [4]:
fname = "CH82.scn" # binary SCN file containing simulated idealised single-channel open/shut intervals
tr = 1e-4 # temporal resolution to be imposed to the record
tc = 4e-3 # critical time interval to cut the record into bursts
conc = 100e-9 # agonist concentration
Initialise Single-Channel Record from dcpyps. Note that SCRecord takes a list of file names; several SCN files from the same patch can be loaded.
In [5]:
# Initaialise SCRecord instance.
rec = dataset.SCRecord([fname], conc, tres=tr, tcrit=tc)
rec.printout()
Plot dwell-time histograms for inspection. In single-channel analysis field it is common to plot these histograms with x-axis in log scale and y-axis in square-root scale. After such transformation exponential pdf has a bell-shaped form.
In [6]:
fig, ax = plt.subplots(1, 2, figsize=(12,5))
dcplots.xlog_hist_data(ax[0], rec.opint, rec.tres, shut=False)
dcplots.xlog_hist_data(ax[1], rec.shint, rec.tres)
fig.tight_layout()
In [7]:
mec = samples.CH82()
mec.printout()
In [8]:
# PREPARE RATE CONSTANTS.
# Fixed rates
mec.Rates[7].fixed = True
# Constrained rates
mec.Rates[5].is_constrained = True
mec.Rates[5].constrain_func = mechanism.constrain_rate_multiple
mec.Rates[5].constrain_args = [4, 2]
mec.Rates[6].is_constrained = True
mec.Rates[6].constrain_func = mechanism.constrain_rate_multiple
mec.Rates[6].constrain_args = [8, 2]
# Rates constrained by microscopic reversibility
mec.set_mr(True, 9, 0)
# Update rates
mec.update_constrains()
In [9]:
#Propose initial guesses different from recorded ones
#initial_guesses = [100, 3000, 10000, 100, 1000, 1000, 1e+7, 5e+7, 6e+7, 10]
initial_guesses = mec.unit_rates()
mec.set_rateconstants(initial_guesses)
mec.update_constrains()
mec.printout()
In [10]:
# Extract free parameters
theta = mec.theta()
print ('\ntheta=', theta)
In [11]:
def dcprogslik(x, lik, m, c):
m.theta_unsqueeze(np.exp(x))
l = 0
for i in range(len(c)):
m.set_eff('c', c[i])
l += lik[i](m.Q)
return -l * math.log(10)
In [12]:
# Import HJCFIT likelihood function
from dcprogs.likelihood import Log10Likelihood
# Get bursts from the record
bursts = rec.bursts.intervals()
# Initiate likelihood function with bursts, number of open states,
# temporal resolution and critical time interval
likelihood = Log10Likelihood(bursts, mec.kA, tr, tc)
In [13]:
lik = dcprogslik(np.log(theta), [likelihood], mec, [conc])
print ("\nInitial likelihood = {0:.6f}".format(-lik))
In [14]:
from scipy.optimize import minimize
print ("\nScyPy.minimize (Nelder-Mead) Fitting started: " +
"%4d/%02d/%02d %02d:%02d:%02d"%time.localtime()[0:6])
start = time.clock()
start_wall = time.time()
result = minimize(dcprogslik, np.log(theta), args=([likelihood], mec, [conc]),
method='Nelder-Mead')
t3 = time.clock() - start
t3_wall = time.time() - start_wall
print ("\nScyPy.minimize (Nelder-Mead) Fitting finished: " +
"%4d/%02d/%02d %02d:%02d:%02d"%time.localtime()[0:6])
print ('\nCPU time in ScyPy.minimize (Nelder-Mead)=', t3)
print ('Wall clock time in ScyPy.minimize (Nelder-Mead)=', t3_wall)
print ('\nResult ==========================================\n', result)
In [15]:
print ("\nFinal likelihood = {0:.16f}".format(-result.fun))
mec.theta_unsqueeze(np.exp(result.x))
print ("\nFinal rate constants:")
mec.printout()
In [16]:
from dcprogs.likelihood import QMatrix
from dcprogs.likelihood import missed_events_pdf, ideal_pdf, IdealG, eig
qmatrix = QMatrix(mec.Q, 2)
idealG = IdealG(qmatrix)
Note that to properly overlay ideal and missed-event corrected pdfs ideal pdf has to be scaled (need to renormailse to 1 the area under pdf from $\tau_{res}$).
In [17]:
# Scale for ideal pdf
def scalefac(tres, matrix, phiA):
eigs, M = eig(-matrix)
N = inv(M)
k = N.shape[0]
A = np.zeros((k, k, k))
for i in range(k):
A[i] = np.dot(M[:, i].reshape(k, 1), N[i].reshape(1, k))
w = np.zeros(k)
for i in range(k):
w[i] = np.dot(np.dot(np.dot(phiA, A[i]), (-matrix)), np.ones((k, 1)))
return 1 / np.sum((w / eigs) * np.exp(-tres * eigs))
In [18]:
fig, ax = plt.subplots(1, 2, figsize=(12,5))
# Plot apparent open period histogram
ipdf = ideal_pdf(qmatrix, shut=False)
iscale = scalefac(tr, qmatrix.aa, idealG.initial_occupancies)
epdf = missed_events_pdf(qmatrix, tr, nmax=2, shut=False)
dcplots.xlog_hist_HJC_fit(ax[0], rec.tres, rec.opint, epdf, ipdf, iscale, shut=False)
# Plot apparent shut period histogram
ipdf = ideal_pdf(qmatrix, shut=True)
iscale = scalefac(tr, qmatrix.ff, idealG.final_occupancies)
epdf = missed_events_pdf(qmatrix, tr, nmax=2, shut=True)
dcplots.xlog_hist_HJC_fit(ax[1], rec.tres, rec.shint, epdf, ipdf, iscale, tcrit=rec.tcrit)
fig.tight_layout()
Note that in this record only shut time intervals shorter than critical time ($t_{crit}$) were used to minimise likelihood. Thus, only a part of shut time histrogram (to the left from green line, indicating $t_{crit}$ value, in the above plot) is predicted well by rate constant estimates.