This function computes the forward log-likelihood of data given a model.
See also:
Up: IndexOur Markov model is a graph with colored vertices. A vertex is called a "state," and its color is a nonnegative integer called its "class." States in the same class are indistinguishable (same amp). To describe the vertices, provide the array
int clazz[Ns] = [class of each state]
int Ns = number of states
Each edge is labeled with its transition rate (probability per second). These form the matrix
Q, a Ns x Ns matrix with
\(Q_{a,b}\) = rate from state a to state b
\(Q_{a,a} = - \sum_i Q_{a,i}\) where \(i \neq a\)
Each \(Q_{a,b} = K0_{a,b} * Ligand_{a,b} * e^{K1_{a,b} * Voltage_{a,b}}\). You provide the Ns x Ns matrices
double **K0, **K1 of kinetic parameters
int **Ligand, **Voltage index of the ligand or voltage signal influencing each rate, or 0
with the diagonals undefined.
A pair of states is either connected in both directions or neither. To indicate un-connectedness, set
\(K0_{a,b} = K0_{b,a} = 0\)
For constant state entry probabilities, as described in (Qin...1996), provide the array
int P0[Ns]
For equilibrium entry probabilities, as described in (Milescu...2005), provide
P0 = NULL
The data consist of one or more parallel signals. The first (index 0) is the Markovian one to be analyzed. Each additional signal describes a ligand or voltage variable. Signals are in model order; i.e., if (v = Voltage[a][b]) != 0, then the signal at index v is the voltage controlling a->b. For each signal you provide
int DwellCount number of events int Classes[DwellCount] class index of each event float Durations[DwellCount] length of each dwell in milliseconds double Amps[ClassCount] amplitude (or ligand/voltage value) of each class
and you give the signals together as
int Nsignal number of signals int *dwellCounts int **classeses float **durationses double **ampses
You provide one or more segments of data, each with the same number of signals. All together they are given as
int Nseg int Nsignal int **dwellCountses [segment][signal] int ***classeseses float ***durationseses double ***ampseses
We multiplex the signals to create a single idealized signal which changes class whenever any source signal changes class. Each plexi-class denotes a Markov class with a specific set of experimental constants (the other signals' idealized amplitudes).
Then, as in MIL, we subtract tdead. tdead is the longest duration of events that can't be reliably detected. MSL deletes any such events, by merging them with their prior. Then for computational reasons, tdead is subtracted from each event, and they're converted to seconds.Sadly, someone has to allocate memory to store this processed signal. It has the form
int Nseg int Nsignal int dwellCounts[Nseg] int **classeses float **durationses int Nplex int plexiclasses[2*Nplex] alternating markov-class[i/2], stimclass[i/2] int Nstim double stimclasses[Nstim*Nsignal] [stimcls * Nsignal + signal_ix], signal 0 undefined
First, we can scan your data to compute upper bounds for dwellCounts, Nplex and Nstim. You allocate dwellCounts[Nseg] and call:
end_html */ extern "C" MAXILL_API void __stdcall msl_multiplex_bounds(int Nseg, int Nsignal, int **dwellCountses, int ***classeseses, float ***durationseses, double ***ampseses, int *out_dwellCounts, int *out_Nplex, int *out_Nstim); /* begin_htmlThen you allocate classeses[Nseg][dwellCounts[i]] and durationses, plexiclasses and stimclasses, and call this function to multiplex and process the signals:
end_html */ extern "C" MAXILL_API void __stdcall msl_multiplex(int Nseg, int Nsignal, int **dwellCountses, int ***classeseses, float ***durationseses, double ***ampseses, double tdead_ms, int *out_dwellCounts, int **out_classeses, float **out_durationses, int *out_Nplex, int *out_plexiclasses, int *out_Nstim, double *out_stimclasses); /* begin_htmlNow you can call the subi_ll function below.
double *LL: upon return, contains log(prob. that model generated this data)
Returns: 0 on success, error codes to be defined.