Waiting time distributions for clusters of complex molecules

Waiting time distributions are in the core of theories for a large variety of subjects ranging from the analysis of patch clamp records to stochastic excitable systems. Here, we present a novel exact method for the calculation of waiting time distributions for state transitions of complex molecules with independent subunit dynamics. The absorbing state is a specific set of subunit states, i.e. is defined on the molecule level. Consequently, we formulate the problem as a random walk in the molecule state space. The subunits can posses an arbitrary number of states and any topology of transitions between them. The method circumvents problems arising from combinatorial explosion due to subunit coupling and requires solutions of the subunit master equation only.

The defining property of self-amplifying reactions and reaction cascades is that a single initial event causes many subsequent reactions and product molecules. When the number of molecules causing the first event is small, initiation occurs stochastically. This is the case in the tiny volume of living cells, where such randomness may even be observable on the cell level [1]. Binding of a ligand to a single G-protein coupled receptor leading to thousands of product molecules downstream in the signalling pathway is such an event [2]. Similarly, the crossing of the activation threshold of an excitable system causes strong amplification of the signal. That crossing may be preceded by a complex process consisting of many random elemental steps. Well-known examples are puffs and sparks in intracellular Ca 2+ release and randomly spiking neurons [1,3,4]. These random processes can be characterised by waiting time distributions, which are inter-puff and inter-spike interval distributions in biological terms.
The above examples have in common that only a few specific events out of the complex dynamics on one structural level drive the dynamics on the next higher level. For instance, among the many state transitions an ion channel may have, only the opening and closing events are relevant for membrane voltage dynamics. This suggests reducing the description of the micro-or mesoscopic level to the states relevant for the dynamics on the macroscopic level and the corresponding waiting time distributions. The complex processes on the microscopic level are subsumed into the distributions. An excitable system defined that way is presented in [5,6]. Thus, waiting time distributions can serve as constituting functions of a theory of hierarchical stochastic systems.
Waiting time distributions have also been applied to model the behaviour of linear molecular motors such as myosin V and kinesin, or of rotary motors like the F1-ATPase [7][8][9][10]. A single motor step is perceived as a sequence of state transitions of the motor molecule. Since these motors always return to the same initial configuration after a stroke, the corresponding waiting time distribution defines a renewal process.
A revived interest in the description of oscillatory systems as a set of discrete states with cyclic, unidirectional state transitions arose recently spurred by the search for universality in the non-equilibrium phase transitions to synchronous oscillations [11] and the quest for the origin of oscillatory dynamics from random elements [5,11]. Waiting time distributions are one of the defining elements of these systems.
In summary, there is a growing interest for characterising systems that arise from direct applications and fundamental questions of physics, in terms of discrete states and waiting time distributions. In this letter, we will calculate waiting time distributions for state transitions of complex molecules which consist of several identical and independent subunits. The activated molecule state is a specific configuration of subunit states. This molecular state is the 38003-p1 absorbing state of our waiting time calculations. Let us assume for the moment that the activated molecule state is the one where all subunits are in the same state a. We consider the event when for the first time all but one of the subunits are already in the state a and the remaining subunit makes a transition into the state a. All individual subunits may have visited the state a many times before all of them meet there. Consequently and although all subunits are independent, the waiting time distribution we would like to compute is not a power of the single subunit waiting time distribution for a transition to a. The absorbing molecule state couples the trajectories of the individual subunits.
Methods based on Fokker-Planck-equations cannot be applied to these problems for two reasons. Firstly, there is no large parameter as, e.g., the volume that allows for a proper derivation from microscopic dynamics, and secondly, the large number of dimensions of the phase space would render a Fokker-Planck equation intractable. We therefore need to solve the master equation for the complete molecule. However, this is often difficult to even write down due to combinatorial explosion of the number of states. Our approach provides a way out of that dilemma since it necessitates to solve the subunit master equations only. We will illustrate our ideas for clusters of molecules, but we start with a single molecule.
This work was inspired by problems from intracellular Ca 2+ dynamics. Calcium release from storage compartments is controlled by ion channels, which we will use to exemplify our concept. The open probability of the channel depends on the calcium concentration on the outside of the storage compartment (see [1] for a review). It increases with increasing Ca 2+ concentration, i.e. Ca 2+ release is self-amplifying. Very high Ca 2+ concentrations inhibit the channel and terminate release. Activation and inhibition together create a bell-shaped dependence of the open probability on [Ca 2+ ].
Release channels are grouped into clusters, each containing 1-40 channels [1]. The spontaneous Ca 2+ release through channels of a single cluster, called a puff, is the elemental random event of intracellular Ca 2+ dynamics. The dependence of the open probability on the Ca 2+ concentration entails that channels in a cluster communicate via the local Ca 2+ concentration. However, this feedback is eliminated when the Ca 2+ concentration is constant throughout a cluster so that channels behave independently. As a good approximation, this is true when all channels are closed or at most one channel is open. The absorbing state can be defined in a way that the waiting time problem describes one of the transitions between these two situations [12].
We consider a channel molecule consisting of h identical subunits. They possess an arbitrary number of states and any topology of transitions between them. The states of the channel molecule can be described as points in a multidimensional state space. Let q(X, t|Y, τ ) denote the probability of the channel to be in state X at time t when it was in state Y at time τ , and F (X, t|Y, 0)dt the probability of arriving in state X for the first time in the interval [t, t + dt] when the initial state was Y . These two quantities satisfy the relation [13] q(X, t|Y, The probability of being in X at time t equals the probability of arriving there for the first time in some interval [τ, τ + dτ ], τ t, times the probability of being in X at time t given the initial condition (X, τ ). The integral includes all τ between 0 and t, and the delta-functions take care of initial values. In the case of a time-homogeneous Markov process, q(X, t|Y, τ ) = q(X, t − τ |Y, 0), and the Laplace transform of eq. (1) for X = Y reads aŝ Since F (X, t|X, 0) = δ(t), we obtain finallŷ The state of a subunit is defined by the binding of ligands to binding sites on the subunit. We model the stochastic association and dissociation of ligands by a master equation. Let n denote the number of states of a single subunit, and p (i, t|j, 0) the probability of it being in state i at time t if the initial state was j. The time evolution of P (t|j, 0) := (p(1, t|j, 0), . . . , p(n, t|j, 0)) is governed by the master equation d t P (t|j, 0) = W P (t|j, 0). The entries w kl , k = l of the matrix W ∈ R n×n represent the transition probability densities from state l to state k. The diagonal elements w ii are such that each column of W sums to 0. We will use 3-state and 4-state subunit schemes in this report. Four states correspond to 2 binding sites on a subunit. Three states could result, e.g., from lumping of the two inhibited states [12] leading to unique rest, activated and inhibited states of a subunit. The transition matrices are given by for the 4-state model and by for the 3-state model. Since our approach holds for any subunit dynamics with constant transition rates, the entries in eq. (4) are chosen arbitrarily. The only restriction is that they comply with detailed balance. Note that the entries have units of inverse time. We denote the 4 states as 00, 10, 11 and 01, where each digit indicates the state of a binding site as occupied (1) or empty (0). The 3 states read as i, a and r. The state of a channel molecule will be described by the numbers of subunits in each of the subunit states, i.e. (n i , n a , n r ) for a channel consisting of 3-state subunits and (n 00 , n 10 , n 11 , n 01 ) for 4-state subunits. With a constant matrix W which can be diagonalised, a solution of the master equation is given by The coefficients c ji ∈ R are determined by initial conditions. λ i ∈ R and V i = (v 1i , . . . , v ni ) denote an eigenvalue and the corresponding eigenvector of W , respectively. In a first step, we focus on channels that are activated when all h subunits are in a unique activated state labelled by the index a. Hence, the activated channel state A is n a = h, n i = 0, i = a. We consider an initial condition I that has m i , i = 1, . . . , n subunits in state i. The probability q := q(A, t|I, 0) of a channel being activated at time t given I at time 0 is kij e lj λj t .
The prime indicates the restriction n j=1 k ji = l i . The ansatz for q in eq. (6) reflects the fact that we consider the subunits to be identical and to gate independently [14]. Since there are r := ( h+n−1 n−1 ) terms in eq. (7) corresponding to r ways to decompose h into n summands -note at the same time that there are exactly r channel states-, we cast eq. (7) into the form q = r j=1 M j (m) exp (η j t) .
η j = l 1 λ 1 + . . . l n λ n is a linear combination of the eigenvalues of W , and M j (m) equals M ({l}, m) evaluated at the j-th decomposition of h. The Laplace transform of eq. (9) reads asq The common denominator ofq(A, s|I, 0) is independent of the initial condition. Therefore, it cancels when we insert eq. (10) into eq. (3). Setting q(A, s|I, 0) : for I = A. The first passage time distribution in the time domain follows readily as However, experiments have demonstrated that channels are conducting in more than one configuration. That is true for the IP 3 R channel in particular, which is deemed open when 3 out of 4 subunits are in the active state [1]. Therefore, we consider the general case that the molecule is activated whenever it is in one of the states A i , i = 1, . . . , r a < r. Let F (A i A, t|I, 0)dt denote the joint probability of being in the state A i and in the activated state A := i A i for the first time in the interval [t, t + dt] given an initial state I at time 0. The probability F (A i A, t|I, 0)dt subsumes all trajectories in phase space that arrive at A i for the first time without having visited any of the other states A j , j = i before. Consequently, the probability density F (A, t|I, 0) for the molecule of arriving in A for the first time is F (A, t|I, 0) = i F (A i A, t|I, 0), since the A i are mutually exclusive [15]. This leads directly to a generalisation of equation (1) as 38003-p3 i = 1, . . . , r a , which is equivalent tô for I = A i . Equation (14) is a system of r a equations for the r a unknownsF (A j A, s|I, 0). Let B ∈ C ra×ra be a matrix with the entries b ij :=q(A i , s|A j , 0), thenF (A j A, s|I, 0) = i B −1 jiq (A i , s|I, 0). To illustrate the preceding scheme, we analyse a molecule that consists of three subunits with three states each. We choose the initial state with all subunits in the third state, I = (0, 0, 3), and the activated state as A = A 1 ∪ A 2 = (1, 2, 0) ∪ (0, 2, 1). This leads to so that, e.g.,F (A 1 A, s|I, 0) = N 1 / det B with The computation ofq(A i , s|I, 0) andq(A i , s|A j , 0) requires an extension of our notion of eq. (6). Now, the transition from A j to A i comprises all trajectories from any subunit configuration in A j to any configuration in A i . Take, e.g., the transition from A 1 to A 2 . Since the subunit that is initially in the inhibited (first) state can either move to the active (second) or to the resting (third) state, we need to consider both possible transitions with appropriate combinatorial factors: q(A 2 , t|A 1 , 0) = p(2, t|2, 0) 2 p(3, t|1, 0) +2p(3, t|2, 0)p(2, t|2, 0)p(2, t|1, 0) . (17) Analogously, q(A i , t|I, 0) must take into account all combinatorial possibilities to get from I to A i , so that, e.g., These state probabilities are sums of exponentials, so that the Laplace transformsq are of the form of eq. (10). It entails thatF (A i A, s|I, 0) are fractions of polynomials of order 2r − 2 and are readily back transformed in analogy to eq. (12). The left panel in fig. 1 depicts the outcome of this concept. The analytical solution is almost indistinguishable from results of stochastic simulations. We now move one structural level up to channel clusters, in order to describe events like puffs in intracellular Ca 2+ dynamics. A puff starts when the first channel in the cluster opens. The probability for the first channel activation depends on the initial states of all channels. Let u := (u 1 , . . . , u r ) with r i=1 u i = N denote the primary distribution of the N channels among the r channel states. Then, assuming that channels in a cluster behave  Fig. 1: Left: Escape probability F for a channel with 3 3-state subunits using eqs. (12) and (14). The initial state is I = (0, 0, 3), and the absorbing state is A = (0, 2, 1) ∪ (1, 2, 0). Right: Escape probabilityF b for a cluster of 7 trimeric channels. Each subunit is modelled by 4 states. 10 is the activated state of the subunit. Initial conditions are (n00, n10, n11, n01) = (2,8,4,7), and the absorbing state is when a channel has 3 subunits in 10 for the first time. The inset shows a blow-up of the delta-peak at t = 0. Lines: analytic results, shaded areas: histograms from stochastic simulations (500000 trials).
independently, the probability distribution for a cluster to be activated for the first time at time t reads as Here, G i is the probability that a channel that was originally in the channel state I i has not yet reached A at time t; note that G(A, t|A, 0) = δ(t). Equation (19) states that a cluster is activated for the first time when one channel opens for the first time, while all other channels are still closed. Since subunits represent the defining dynamical elements, initial data are often given as the number of subunits in each of the n states. However, this does not determine the distribution of the Nh subunits among the N channels. Consequently, we need to average over all possible initial cluster states with b i , i = 1, . . . , n subunits in the i-th state. This leads to the probability distribution F b (t) = (1/Γ) u f u F u (t), where the sum runs over all r-tuples u obeying and where The constants m ij represent the number of subunits in the state i given the j-th channel state ( i m ij = h). Hence, a single square bracket equals the number of possibilities to generate a distinct channel state, so that f u corresponds to the number of ways to realise a cluster state where u i channels are in the state i. The normalization, Γ := u f u , can be shown to evaluate to (Nh)!/ (b 1 ! · · · b n !). The right panel in fig. 1 compares results forF b to stochastic simulations of an entire cluster. The agreement between analytic results and simulations is excellent. Note the resolution of the delta peak at t = 0. Although the graph illustrates only small fractions of the probability distributions, the concord persists throughout the tails.
In summary, we have presented an exact method to calculate waiting time distributions for clusters of stochastic molecules from master equations. Molecule subunits should be coupled by the absorbing state only but do not need to be identical (as we assumed for convenience here) in order for the method to be applicable. Constructing the Laplace transform of the waiting time distributions from the Laplace transform of solutions of the master equation of single subunits circumvents combinatorial explosion arising from the coupling of the subunit master equations by the absorbing state. Note, that the method even calculates splitting probabilities for first passage to different absorbing states along the way.
New approaches for solving non-factorisable master equations have been suggested in recent years. They mainly focus on numerical schemes and approximations or involve complex computation for subunits with more than two states (see, e.g., [16,17]). The advantage of our method lies in utilising explicitly the independence of the subunits despite the coupling by the absorbing state. Hence, we start directly from the master equation of a subunit, and not from the full molecule master equation. This drastically reduces the dimension of the phase space.
An effective method to calculate waiting time distributions, as we have demonstrated here, paves the way to a description of hierarchical stochastic systems. Such a theory is needed in particular for the modelling of processes in living cells, where we have small numbers of active elements that exhibit discrete states.