July 26, 2008

The Q Continuum

The week's course turned out to be rather interesting, although the kind of analysis under discussion is somewhat more detailed than I'm likely to need anytime soon. One thing that soon became clear was that my own level of engagement was inversely proportional to how software-related things happened to be. The underlying mathematics, which appears so elegant in retrospect but must have been an absolute bastard to derive, always kept me gripped. DC's slightly rambling, donnish explication of it was a privilege to experience. The implementation, whether pared-down in MathCad for educational purposes or incarnated in what apparently represents the scientific state of the art in DC's program suite... um, not so much.

The thing is, I am almost entirely persuaded that the techniques under discussion represent the right -- which is to say, the best available and most theoretically correct -- way to go about things. I can also see why so few people do. It's not because the whole thing is conceptually difficult -- bits of it certainly are, but not outlandishly so in the already-abstruse context of electrophysiology. But it simply isn't available in a form that is accessible to someone not already institutionalised by years of exposure to the software quirks, even if that someone is otherwise well-versed in the intricacies of single channel recording.

At one level, I'd really like to do something to improve this particular situation, but it's not clear what that might be. The obvious approaches represent a huge amount of work, and would be to the benefit of almost no-one; they are, in consequence, not all that appealing. But I'll be looking into it in the coming weeks, so maybe there'll be some worthwhile contribution I can make.

I have, you may note, completely glossed over the actual meat of the subject, and I don't propose to attempt any kind of explanation here; at least, not just now. But, purely as a token, here's the "mug" equation, which gives the expected burst length distribution for an ion channel whose behaviour is characterised by a mechanism Q, assuming the recording apparatus is perfect, which it never is:

fb(t) = φb [eQEEt]AA (GABGBC + GAC) uC

And here's some Matlab code to calculate this, along with the corresponding open and shut time distributions, under the same unrealistic measurement conditions.

function [A lambda] = spec ( M )
% SPEC generates the spectral expansion matrices of the supplied matrix
% [A lambda] = SPEC(MATRIX)
% A is a cell array of the spectral expansion matrices
% lambda is a vector of the corresponding eigenvalues
[h v] = size(M);
if ( h ~= v ), error 'supplied matrix must be square!'; end

A = {};
[V, D] = eig(M);
rV = inv(V);
lambda = diag(D);

for k = 1:h
    cE = V(:, k);
    rE = rV(k, :);
    A{k} = cE * rE;

function B = sylv ( M, f, args )
% SYLV calculates a function of a matrix using Sylvester's Theorem
% FUNC must be of the form RESULT = F(MATRIX, ARG)
% B is a cell array containing the results for all ARG in ARGS
[h v] = size(M);
if ( h ~= v ), error 'supplied matrix must be square!'; end

[A lam] = spec(M);
B = {};

for r = 1:length(args)
    Br = f(lam(1), args(r)) * A{1};
    for k = 2:h
        Br = Br + f(lam(k), args(r)) * A{k};
    B{r} = Br;

function [D_o D_s D_b] = qideal ( Q, kA, kB, kC, t )
% QIDEAL calculates the ideal open, shut and burst time distributions for a Q matrix
% [f_open f_shut f_burst] = QIDEAL( Q, kA, kB, kC, t )
% Q must be in canonical form, with A (open) states before B (fast shut) before C (slow shut)
% kA, kB and kC are the number of each kind of state
% t is a vector of time points for which to calculate the f values
% no account is taken of the distortion caused by missed events, and error checking is minimal
[h v] = size(Q);
if ( h ~= v ), error 'Q matrix must be square'; end

kF = kB + kC;
kE = kA + kB;
k = kA + kF;

if ( k ~= h ), error 'kA + kB + kC is not the same as k'; end

% get the partition matrices
Q_AA = Q(1:kA, 1:kA);
Q_AF = Q(1:kA, (kA+1):k);
Q_FA = Q((kA+1):k, 1:kA);
Q_FF = Q((kA+1):k, (kA+1):k);

Q_BB = Q((kA+1):(kA+kB), (kA+1):(kA+kB));
Q_CB = Q((kE+1):k, (kA+1):kE);
Q_CA = Q((kE+1):k, 1:kA);
Q_BA = Q((kA+1):kE, 1:kA);
Q_AB = Q(1:kA, (kA+1):kE);
Q_EE = Q(1:kE, 1:kE);

% define the (time independent) G matrices
G_AF = inv(-Q_AA) * Q_AF;
G_AB = inv(-Q_AA) * Q_AB;
G_BA = inv(-Q_BB) * Q_BA;

% calculate equilibrium occupancies
S = [Q, ones(k,1)];
p_inf = ones(1,k) * inv(S * S');

% calculate initial state vectors
pF_inf = p_inf((kA+1):k);
phi_o = pF_inf * Q_FA;
phi_o = phi_o / sum(phi_o);

phi_s = phi_o * G_AF;

pC_inf = p_inf((kE+1):k);
phi_b = pC_inf * (Q_CB * G_BA + Q_CA);
phi_b = phi_b / sum(phi_b);

% calculate burst-end vector
e_b = (eye(kA) - G_AB * G_BA) * ones(kA, 1);

% f_open(t) = sum(phi_o * exp(Q_AA * t) * Q_AF)
% f_shut(t) = sum(phi_s * exp(Q_FF * t) * Q_FA)
% f_b(t) = phi_b * exp(Q_EE * t)_AA * (-Q_AA) * e_b
f = @(x,t) exp(x*t);
expQA = sylv(Q_AA, f, t);
expQF = sylv(Q_FF, f, t);
expQE = sylv(Q_EE, f, t);

D_o = zeros(size(t));
D_s = zeros(size(t));
D_b = zeros(size(t));

for ti = 1:length(t)
    D_o(ti) = sum(phi_o * expQA{ti} * Q_AF);
    D_s(ti) = sum(phi_s * expQF{ti} * Q_FA);
    tmp = expQE{ti};
    D_b(ti) = phi_b * tmp(1:kA,1:kA) * (-Q_AA) * e_b;

Generating realistic distributions is a much more difficult problem. If I manage to get to grips with the key 1990 and 1992 papers by Hawkes, Jalali & Colquhoun I may post code for that too. I want it clear from the start, however, that I will not attempt any kind of fitting. That stuff has an SEP field the size of Kazakhstan around it...
Posted by matt at July 26, 2008 10:16 PM

Something to say? Click here.