Single molecule dynamics in a virtual cell: a three-dimensional model that produces simulated fluorescence video-imaging data

The analysis of single molecule imaging experiments is complicated by the stochastic nature of single molecule events, by instrument noise and by the limited information which can be gathered about any individual molecule observed. Consequently, it is important to cross check experimental results using a model simulating single molecule dynamics (e.g. movements and binding events) in a virtual cell-like environment. The output of such a model should match the real data format allowing researchers to compare simulated results with the real experiments. The proposed model exploits the advantages of ‘object-oriented’ computing. First of all, the ability to create and manipulate a number of classes, each containing an arbitrary number of single molecule objects. These classes may include objects moving within the ‘cytoplasm’; objects moving at the ‘plasma membrane’; and static objects located inside the ‘body’. The objects of a given class can interact with each other and/or with the objects of other classes according to their physical and chemical properties. Each model run generates a sequence of images, each containing summed images of all fluorescent objects emitting light under given illumination conditions with realistic levels of noise and emission fluctuations. The model accurately reproduces reported single molecule experiments and predicts the outcome of future experiments.


Introduction
The past two decades have been marked by the rapid progress in single molecule imaging in live cells. In the early studies, small particles, specifically attached to membrane molecules, were used to track movements of individual molecules [1][2][3]. Specific fluorescent probes, containing tens of molecules, were also used for similar purposes [4,5]. Progress in camera technology, optics and lasers made it possible to detect and track single membrane-associated molecules tagged with a single fluorophore [6,7] or fused to a single fluorescent protein [8,9] which could be reliably identified as the single molecule of interest. Fast sensitive cameras have allowed detection and tacking of large protein complexes moving within the cytoplasm and at the cell nucleus [10,11]. The application of total internal reflection fluorescence microscopy (TIRFM) has allowed researchers to illuminate a very thin layer of the solution close to the surface of the glass-water interface which greatly improves the signal-to-noise ratio [12], critical for single molecule detection. Some intracellular proteins containing membrane binding domains were detected and tracked at the plasma membrane, because molecules slow down and become visible as discrete spots of light upon binding to the membrane [6,[13][14][15]. Some fine details of molecule movements, which include confined diffusion [2,16], presence of membrane barriers [17] or membrane microdomains [7,18], were evaluated using single molecule or single particle-tracking methods.
One of the most important issues in single molecule research is the stochastic nature of single molecule events, which require large datasets to be analysed in order to obtain statistically valid conclusions. This problem can be solved by using automatic detection and tracking algorithms [5,19] that eliminate operator bias associated with the laborious manual processing of data.
Another important issue is a requirement to cross-validate the results obtained by automated analysis procedures. This problem can be addressed by using computer models to simulate the mobility and binding kinetics of the molecules of interest and to generate the results in the same format as real data. Many of the published works contain a modelling section [10,11,17,[19][20][21] simulating experimental results. Naturally, these models are restricted to a particular problem and are limited in structure and scope. Andrews & Bray [22] proposed a general model emulating chemical interactions between individual molecules of one type moving inside a virtual bacterial cell. This model was modified by Tournier et al. [23], who proposed calculating the probability of interaction between two molecules separated by a significant distance, so that the model could use much longer time steps than the Andrews-Bray model. The models can produce results in the form of spatial distributions of individual objects in a particular chemical state and as a time course of the concentrations of the reacting species. These models use only one type of molecule (free moving in cytoplasm), require a significant number of calculations, and do not simulate imaging conditions affecting the results of real experiments.
The purpose of this work was to construct a novel computer model that can simulate a few distinct 'classes' of single molecules moving both within the cytoplasm and at the cell membrane and simulate chemical interactions between molecules of the same or few different classes. The model takes into account the three-dimensional illumination pattern created under given conditions (epi-illumination, TIRFM, confocal microscopy) which would affect the emission rate of individual molecules. The images of all the lightemitting molecules, within and beyond the focal plane, are projected onto a virtual imaging device (e.g. EMCCD camera) using rules of optics and experimentally or empirically determined noise and signal characteristics. The model is optimized to use the minimal number of calculations and has a modular structure giving it the ability to be extended to more complex scenarios such as abnormal diffusion, directed movements and single molecule dynamics in the presence of some intracellular structures (e.g. nucleus or cytoskeleton), and to add any new properties (e.g. take into account the effects of polarized illumination and the orientation of the molecules).
The proposed model was extensively tested under a number of scenarios. The results of the modelling (analysed with the use of the automatic detection and tracking algorithms [19]) accurately reproduce the results of the reported experiments [13,14,[24][25][26][27][28]. The model was also used to evaluate the putative effects of partial permeability of membrane barriers and of the viscosity of lipid rafts on the outcome of single molecule experiments.

Model
The basic model may contain objects of one or two major classes: cytoplasm-based molecules moving inside an enclosed volume ('cell'), and/or objects moving on the surface of this volume ('plasma membrane'). An object of either class has a set of default properties including physical coordinates (x, y, z); mobility coefficient; indexes for a bound object(s) of the same or other classes; up to eight fluorescent tags; a colour of fluorescent tag (if required). Intracellular 'molecules' can move freely inside a 'cell', and (if permitted) bind to the 'molecules' belonging to 'membrane-based' classes. The 'bound-tomembrane' intracellular molecule moves together with its membrane-embedded 'partner' until the moment of dissociation. It is assumed that the mobility of such a pair is equal to the mobility of the slow-moving membrane-based molecule because it is embedded in a viscous membrane. If two membrane-embedded molecules with mobility D a and D b bind each other, the mobility of the pair (D pair ) will be proportionally slower than the mobility of each member of the pair: D pair ¼ (D a )/(1 þ D a /D b ). For example, if two membrane molecules form a homodimer (see §2 in the electronic supplementary material), then its mobility will be two times slower than the mobility of the monomer (see table S4 in [24]).
During each iteration cycle of the model (figure 1), the physical positions of all the objects are changed according to the motion rules (e.g. random walk, directed movement or static) and the duration of a time step. The possible events, such as binding and dissociation, are executed and the properties of the involved objects are updated. If the 'cell' is illuminated at a given time step, then a fluorescent image is built by 'projecting' the images of all active (nonbleached) fluorescent molecules to the surface of a virtual imaging device. The size and the brightness of a single fluorophore image will depend on its point spread function (PSF) and on its x-, y-, z-position, which determines the amount of incident light the object receives. If required, a number of sequential time steps can be summed to create an averaged fluorescence image which would better represent fastmoving objects. Some random number of normally  Figure 1. Model flowchart. Basic time step includes calculations of all the movements and possible binding -dissociation events for all the objects present in the system. The x-, y-, z-coordinates of all active (light-emitting) objects are used to create a fluorescent image of a cell (if illumination is 'ON'). The number of images can be summed into a single camera frame to create a realistic image of fast-moving molecules.
distributed 'noise counts' can be added to every pixel in the image to simulate background fluorescence and camera noise. The 'raw' fluorescence image can be converted into a synthetic EMCCD image by multiplying the value of every pixel by a random number. The average of these normally distributed random numbers is determined by a camera 'gain'. The synthetic biological 'cell' is represented as an enclosed volume in a shape of a parallelepiped with arbitrary x, y, z sizes (figure 2a). The simulated 'cell' is placed on a 'coverslip' (an ellipse in figure 2a). The interface between the 'cell' and the 'coverslip' defines the focal plane at z ¼ 0.

Object movements
The objects inside the 'cell' normally undergo an unrestricted 'random walk' type of diffusion limited by the 'cell' boundaries. Object movements in all three dimensions (x, y, z), are modelled independently [22] using a Gaussian random number generator (GRNG), based on the Box-Muller algorithm (default output d ¼ 1, m ¼ 0) [29]. Mean-squared displacement (MSD) in one dimension is calculated as 2D D t, where D is the diffusion coefficient, and Dt is a time parallelepiped step [30]. The root-mean-squared displacement (RMSD) was used as a multiplier in the GRNG to generate individual displacement values. The object's x-, y-, z-coordinates were stored as floating-point values to avoid bias caused by rounding of integer numbers.
Objects moving on the cell 'membrane' (another class) undergo two-dimensional random walk. The same algorithm as above was used to calculate the object displacements, but only in two dimensions. Note that each random displacement value was generated separately for each dimension, so if an object, moving on one side of a 'cell', reaches an edge, then it continues its movement on the adjoining side of the parallelepiped (figure 2a). Molecules bound to each other have the same physical coordinates (x, y, z) and move together each time parallelepiped step, so that the displacement values are calculated once for the pair of objects. The changes in the mobility pattern of the paired objects are explained in §2.

Interactions between the objects
Molecules of the same or different species can be allowed to bind each other. Binding events could happen when the corresponding 'binding sites' on both objects are unoccupied and the number of collisions between the two objects reaches a certain threshold. Collisions would occur if two molecules move close to each other, so that the distance between their centres became smaller than the interaction distance (ID) or reaction radii [22]. The simplest way to simulate this process is to reduce the model time step, so that the average amplitude of movements is smaller the ID. Molecules separated by a distance less than or equal to the ID will have at least one collision per time step. If we increase the time step, then the objects will move longer distances, and we will underscore the collision events, because we will not know whether the trajectories of the objects were within the ID or not during the time step.
A practical example. A membrane-bound molecule with D ¼ 0.1 mm 2 s 21 has RMSD % 6.3 nm at Dt ¼ 100 ms, and %2 nm at Dt ¼ 10 ms (MSD ¼ 4D D t). Therefore, a 10 ms time step is sufficient to simulate binding processes at ID ! 5 nm. However, molecules move much faster within the cytoplasm, where protein D is in the range 2 -20 mm 2 s 21 [10,11,31], and, for example, with D ¼ 5 mm 2 s 21 , RMSD % 17.3 nm at Dt ¼ 10 ms, and %5.5 nm at Dt ¼ 1 ms Thus, the above-mentioned approach [22] will consume substantial computational time, because we might want to model collisions of hundreds or even thousands of molecules over a reasonable period of time (e.g. 10 -100 s). Another approach, proposed by Tournier et al. [23], uses probabilitybased description of the reaction process, suitable for the longer time steps. A simplified version of this approach, requiring a minimal number of calculations, was developed for the present model.
To calculate the probability of binding between any potential pair of objects we need to calculate the average interaction time (AIT)-that is, the time which this couple would spend within the ID during the current time step. For example, we can assume that molecules A and B are point-like objects randomly moving in one-dimensional space; A and B are separated by a distance x; we can consider molecule A as a fixed target, whereas molecule B has a mobility of If we assume that d x ¼ ID we can calculate the probability of finding molecule B within the ID of molecule A at time t. We can use discrete integration of the above function over time to find the AIT for any given x. Similar equations can be used to calculate the AIT for objects moving in two-and three-dimensional space, but, in these cases, we replace the ID with interaction area (IA ¼ pr 2 ), and interaction volume and Integrating the probability density function for all possible pairs of objects each time step is unnecessary. Instead, we can build look-up tables containing the AIT values for each class of molecule allowed to bind each other. The index of the look-up table (integer value in 'nm' units) can be used to match the distance x separating objects at the beginning of a time step and AIT. These tables need to be recalculated only when we change the mobility, or ID, or a time step. The example in figure 3a shows how AIT depends on mobility of the membrane-bound objects (twodimensional diffusion) at a time step of 100 ms (ID ¼ 5 nm). Figure 3b shows how the AIT would depend on the duration of a time step if, for example, we simulate dimerization of the membrane-bound molecules. It is clear that the AIT decreases dramatically when the initial distance between molecules increases or when we reduce a time step. In many cases, it will be unnecessary to check the probability of binding for the objects placed at a distance greater than 1 mm apart. Examples in figure 3c,d show the AIT for the molecules moving in a three-dimensional space.
Because we know the AIT, we can use a 'binding rate coefficient' or 'reaction rate' R bind (s 21 ) to calculate the probability of binding P bind ¼ AIT . R bind . It was assumed that R bind is always smaller than the diffusion-limited binding rate [22,32,33], especially in the case of ligand binding where it depends on the orientation of the small binding pocket at the moment of collision [33], so that more than one collision is required for a successful binding. Therefore, we can use a linear random number generator (with output proportional to P bind ) to calculate the outcome of a binding test for any potential pair of molecules.
Bound molecules dissociate according to a simple zeroorder chemical kinetics [32], which is modelled as a random process. The probability of dissociation was calculated as The R bind and R diss can be changed during a model run to simulate, for instance, some real responses to biological stimulation.

Illumination of a specimen
In order to calculate the number of photons (N ph ) emitted by a fluorescent tag during each time step and calculate the probability of bleaching we need to know the illumination intensity at the x, y, z location of an object. There are three major illumination modes (figure 1e) which are used in the present model: epi-illumination, TIRFM illumination and confocal illumination.

Laser epi-fluorescent illumination
Laser epi-fluorescent illumination usually uses a TEM 00 laser beam which has a two-dimensional Gaussian profile (figure 2e, left panel) measured by its full width at half maximum (FWHM). It is assumed that the illumination intensity I(x,y,z) is constant in a z-direction, because we are dealing with a transparent specimen. We can use a two-dimensional Gaussian function to calculate intensity I(x,y,z) in any point of a specimen, where d is a standard deviation (for Gaussian function

TIRFM illumination
TIRFM illumination exploits the effect of total internal reflection of light occurring at the interface between high and low refractive index materials when the angle of incident light u is greater than 'the critical angle' u c ¼ sin 21 (n 1 /n 3 ), where n 1 and n 3 are the refractive indices of the low-and high-index materials, respectively [12].
If we assume that a TEM 00 laser beam is used for TIRFM illumination, then it will have a two-dimensional Gaussian profile in the x,y-plane (see epi-illumination mode described above), but the illumination intensity will decrease exponentially in the z-direction (figure 2e, middle panel) above the interface between the two media. The intensity profile in the z-direction will depend on the angle and wavelength of the illuminating beam and the refractive indices of both media [12], where l 0 is the wavelength of the incident light in a vacuum, u is the angle of incidental light, and n 1 and n 3 are refractive indexes of the liquid and the solid, respectively. Depth d is independent of polarization of an incident beam and decreases with increasing u. Except for u ! u c (where d ! 1), d is of the order of l 0 or smaller [12]. We can calculate the intensity at any point of a specimen under TIRFM illumination as where I (0,0,0) is an intensity in the centre of the field at the coverslip level, and d is a standard deviation. The flatness of the x,y-plane of illumination in epiillumination and TIRFM modes depends on the FWHM of a TEM 00 laser beam. The profile can be virtually flat if the beam FWHM ) the x, y sizes of the illuminated area, but if the laser beam is narrow (FWHM the x, y sizes of a cell) then the illumination profile will be uneven, which can severely affect the results of measurements (see example in figure 4a). Note that if the simulated laser beam differs from the TEM 00 mode beam, other functions should be used to describe its intensity profile in the x, y plane. The diameter of the illumination field can be 'truncated' by a circular 'diaphragm' or by a mask of any other shape or size.

. Confocal illumination
Confocal illumination ( figure 2e, right panel) is the simplest to model because it is a scanning method, where every effort is made to ensure that illumination intensity I (x,y) is the same at each x,y-point illuminated during the scan cycle. It is assumed that a high numerical aperture objective lens (NA ! 1.4) is used for imaging, so that the illumination intensity is reduced above and below the focal plane as I(z) ¼ I(z f ) . (Sz f /(Sz f þ p . (tan(u) . |z 2 z f |) 2 ), where I(z f ) is the illumination intensity at focal plane z f , Sz f is the size of the focused beam at z f , and u is half of the angle of the illuminating cone of light (figure 2e, right panel) defined by a NA of the simulated objective lens.

Building the fluorescent image
The image of the cell is generated by summing emission from all the fluorescent molecules. In general, the photons emitted by each single fluorophore ( point source of light) are spread according to a realistic PSF. If we assume that our system is equipped with an ideal objective lens, then the PSF of a single fluorophore will be approximated by a two-dimensional Gaussian function (equation (2.4)), where I (0,0) is the intensity in the centre of the spot. FWHM of a point source of light is roughly equal to half of a wavelength of the simulated colour (FWHM ¼ 0.6l/NA, where l is a wavelength and NA is the numerical aperture of the objective lens). The average number of photons (N ph ) emitted by a single fluorophore is determined by the illumination intensity at its location: N ph(x,y,z) ¼ I xyz . N ph(0,0,0) . A linear random number generator is used to simulate the 'photon noise' characteristic for single fluorophore emission, so that the range of N ph fluctuations is (N ph ) 1/2 . If an object moves above or below the focal plane of an objective lens, then the size of its image increases. GRNG was used to generate x-and y-coordinates of each of the N ph photons accumulated in the corresponding 'pixels' of a virtual CCD camera. The FWHM value used as the GRNG multiplier was determined by the distance between the z-coordinate of an object and a focal plane. In the case of an ideal objective lens, the FWHM value would increase according to a simple hyperbolic function: FWHM(z) ¼ g obj . (z 2 z f ) 2 þ FWHM f , where g obj is a coefficient defined by the magnification of an objective lens and its numerical aperture, z f is the focal plane, and FWHM f is the size of a fluorophore image at z ¼ z f . g obj is equal to 0.6 mm 21 for the objective lens Â 100 NA 1.45 (see §1 in the electronic supplementary material).
In the case of confocal imaging, a significant fraction of photons, emitted by fluorescent molecules placed above or below the focal plane, does not reach the detector because they are blocked by a pinhole. Therefore, these molecules also generate a smaller number of photons than the molecules placed at the focal plane: where N ph (z f ) is an average number of photons detected at the focal plane z f , and g PH is a coefficient determined by the size of the pinhole and the magnification of the objective lens.
In real experiments, some fluorescent molecules will be bleached during the observation period. The photobleaching process is modelled as a simple zero-order chemical reaction, where the probability of bleaching P bleach ¼ R bleach . Dt . I xyz , where R bleach is the photobleaching coefficient (s 21 ) and Dt is the model time step. R bleach can be set before and changed during the model run. Note that if a single molecule object has few fluorescent tags, then each of them would undergo the bleaching process independently of others.

Single particle tracking
The sequences of fluorescent images were analysed by the single particle-tracking software GMIMPRO [19]. The results of automatic detection and tracking were presented as individual trajectories yielding information about mobility and intensity of each individual object. Some single molecule objects were not detected because their trajectories were shorter than the minimal trajectory length (20 data points) set during the detection phase of the analysis. This would happen because the images of individual moving molecules overlapped, or because molecules were bleached, or because molecules moved to the unilluminated parts of the cell (e.g. above the basal side of the cell illuminated in TIRFM mode).

Default conditions
In the examples presented below, the model was run with the cell sizes ¼ 10 Â 10 Â 10 mm 3 , imaging rate 33 fps (10 model time steps summed per frame). An objective lens Â 100/1.45 NA was simulated giving a final magnification of 100 nm pixel 21 . Single fluorophore image size (FWHM) was set to 250 nm, emission rate 6000 photons s 21 in the centre of the illumination area at z ¼ 0, where the photobleaching rate R bleach ¼ 0.5 s 21 .

Random diffusion in cytoplasm and at the plasma membrane
In the simplest case, the model would contain only one class of objects. For example, a single frame from the electronic supplementary material, video S1 (figure 2b) shows randomly moving non-interacting fluorescent molecules (e.g. green fluorescent protein) inside the cell (epi-illumination mode, concentration 2 nM). Electronic supplementary material, video S2 shows the simulation made under the same conditions, but in TIRFM illumination mode. The focus plane was set at z ¼ 0 mm in the first half and then moved to z ¼ 0.2 mm in the second half of the record. Another basic case scenario is molecules randomly moving on the cell membrane without interactions with other molecules or obstacles. Electronic supplementary material, video S3 shows the results of a simulation carried out in epi-illumination mode at a density of 1 molecule mm 22 . In the first half of the record, the focal plane was set to z ¼ 0 mm (figure 2c) and it was then increased to z ¼ 1 mm in the second half of the record (figure 2d). Unlike the previous case, individual membrane-bound molecules were present on the image sequence for some significant time. Electronic supplementary material, video S4 was made under the same conditions, but in TIRFM mode. The imaging conditions can affect both the appearance of fluorescent molecules (brightness and shape) and the results of their detection. Figure 4a shows how the shape of the illuminating beam can affect the results of detection and tracking. When a narrow laser beam was used to illuminate membrane-bound molecules (figure 4a, black line and left inset), the distribution of intensities of tracked molecules was much wider than with the nearly-flat illumination conditions (FWHM ¼ 200 mm). In the last case, the results of automatic tracking had a narrow distribution of intensities of detected and tracked objects (figure 4a, grey line). This conclusion is correct for both TIRFM and epi-illumination microscopy modes. The averaged mobility of automatically tracked objects was 0.298 mm 2 s 21 and the averaged MSD -Dt plot had a linear shape (figure 4b, black line) characteristic of freely diffusing molecules. The grey line in figure 4b shows the averaged MSD-Dt plot for a record made under the same conditions, but the mobility was reduced to 0.01 mm 2 s 21 . The two exponentially decaying lines (right y-axis) in figure  4b show the number of objects which contributed data to the MSD-Dt plot at the increasing time intervals. The fitted off-rate or the rate of losing tracked objects was much higher for the fast-moving objects (1.25 s 21 at D lat ¼ 0.3 mm 2 s 21 ) than for the slow-moving objects (0.55 s 21 at D lat ¼ 0.01 mm 2 s 21 ), where the off-rate was close to the photobleaching rate (R bleach ¼ 0.5 s 21 ). This example shows that even under simple conditions the results of analysis (e.g. calculated bleaching or dissociation rate) can be severely affected by the imaging conditions-overlapping of the images of fast-moving molecules leads to erroneous truncation of the trajectories of detected objects. Choosing cells with lower densities of molecules or decreasing the temperature [28] would allow us to measure the off-rate free from the artefacts of overlapping moving objects.

Anomalous diffusion at the plasma membrane
The movement of molecules associated with the plasma membrane can be affected for a number of reasons. One type of potential restriction is the presence of 'barriers' or 'fences' creating compartments at the plasma membrane [2,16,17,38]. In this case, the slope of the MSD-Dt plot will decrease at higher Dt and become horizontal at MSD values close to 40% of the size of a compartment (see equations 11-13 in [2]). This scenario was simulated by dividing the cell membrane into square compartments (figure 5a, top inset), using virtual 'barriers' restricting the objects' movements into any neighbouring compartment. The example in figure 5a,b shows that if the barriers are completely impenetrable, the average mobility decreases twofold (figure 5a, grey line, and the electronic supplementary material, video S5) compared with the freely diffusing molecules described previously (figure 4b). The MSD-Dt plot became horizontal at MSD values close to 40% of the size of the compartment. However, if the barriers became partially permeable, the MSD-Dt plot became more linear (figure 5b), which makes it difficult to distinguish between the slow diffusion and fast diffusion limited by the barriers.
Another potential factor limiting movement at the cell membrane is the presence of lipid rafts-that is, patches of membrane that have high viscosity compared with the  [13,19,24,28] barriers at cell membrane TIRFM grid sizes 1 Â 1 mm 2 S5 [2] lipid rafts on cell membrane TIRFM raft sizes 1 Â 1 mm 2 , distance between rafts 1 mm  surrounding membrane [7]. The simplest case of a lipid raft scenario is shown in figure 5c, top inset, where square patches are placed at equal distances from each other. On the example shown in figure 5c,d, 1 Â 1 mm 2 patches were separated by 1 mm gaps, so the rafts cover 25% of the cell membrane. The single molecule objects randomly seeded on the membrane at the beginning of the run were concentrated in the raft areas during the run. The rate and degree of crowding would depend on the ratio of the mobility inside/outside the raft patches (R raft ) and the absolute mobility value. If the ratio was large (R raft ¼ 0.1), then we would observe rapid accumulation of molecules in raft patches (figure 5c, middle inset, and the electronic supplementary material, video S6).
After an initial period of equalizing, the distribution of the mobilities of tracked molecules became bimodal (figure 5c, grey line), reflecting the fact that the majority of molecules were concentrated inside the rafts [7]. It is important to note that many molecules inside the crowded patches were not detected and therefore were excluded from the data analysis. The effect of the presence of the rafts was less detectable when the mobility inside the rafts was close to the mobility of the surrounding membrane (R raft ¼ 0.8; figure 5c, black line). The shape of the averaged MSD-Dt plot remained linear (figure 5d) in a wide range of mobility ratios (0.05-0.8). Thus, analysis of the distributions of mobility and detection of clusters of molecules on the membrane is the best method for detecting potential lipid rafts or other types of membrane patches with increased viscosity.

Binding kinetics
Many classes of intracellular molecules can bind to the plasma membrane and become membrane bound for some period of time [10,14,15,31]  specific phosphoinositol phospholipids at the plasma membrane [14]. Because PH molecules bound to a myoblast's membrane have very limited mobility, it was possible to use the gated-illumination time-lapse imaging mode to measure the dissociation rate. The binding rate (R bind ) was set to 1 Â 10 6 s 21 and R diss ¼ 0.05 s 21 . The concentration of fluorescent PH-domain molecules was 2 nM, and the density of nonfluorescent 'phosphoinositol' targets was 1 mm 22 . Under these conditions, about 15% of PH-domain molecules were bound to the plasma membrane (figure 6a and the electronic supplementary material, video S7). To measure the dissociation and photobleaching rates, the model was run in the time-lapse mode (figure 1) with increasing dark intervals between recorded images. The off-rate (the rate of disappearance of the fluorescent spots 'landed' on the membrane during recording) was measured and plotted against the illumination duty ratio. The slope of the fitted regression line was a measure of the photobleaching rate (0.499 s 21 ), and the intercept of the regression line with the y-axis was a measure of the dissociation rate (0.046 s 21 ) [14]. Some molecules, moving freely on the plasma membrane, can bind to other molecules of the same (e.g. form dimers [6,24,35]) or other (e.g. bind to immobile anchors [3,25,27,34]) species. The last case scenario can be modelled by the introduction of the two classes of membrane-bound molecules (one mobile and one immobile) and by allowing molecules of one class to bind molecules of another class. For example, randomly moving membrane-bound fluorescent molecules (density 0.5 mm 22 , D lat ¼ 0.3 um 2 s 21 ) bind to immobile non-fluorescent molecules (density 0.5 mm 22 ) randomly placed at the cell membrane (electronic supplementary material, video S8). The proportion of bound/free moving molecules will depend on the binding and dissociation rates. At R bind ¼ 1 Â 10 5 s 21 and R diss ¼ 1 s 21 , 45% of mobile molecules were bound to immobile counterparts (approx. 10% at R bind ¼ 1 Â 10 4 s 21 ). The results of detection and tracking, at R bind ¼ 1 Â 10 5 s 21 , showed a bimodal distribution of the mobilities of the detected objects (figure 6b), where 23% of objects were completely immobile (average D lat ¼ 0.16 um 2 s 21 ). The fraction of tracked objects which showed a detectable period of immobilization (mobile-immobile-mobile pattern) was small (358 of 4113, or approx. 9% of all objects detected in 10 model runs). The distribution of 'bound' time had an exponential shape (fitted off-rate 1.2 s 21 ; figure 6b, inset). Because after the moment of dissociation the tracked molecules were still visible, we do not have to take into account the effect of photobleaching.

Discussion
This model uses the advantages of object-oriented programming [29,39] to model diverse populations of single molecules represented as objects of a given class. This approach makes it easy to modify and extend the basic model in order to simulate new types of experiment (for example, creating new 'child' classes of single molecules moving on filaments from the 'parent' class of cytoplasm-based molecules) or to use specific parts of the model (e.g. one specific class and one type of imaging condition). The continuous space model, which uses floating-point physical coordinates, was chosen over the matrix-type model [34,37] because it was not limited by the size of the matrix elements in a three-dimensional space and, therefore, allowed modelling of very small changes in the object's position (less than 1 nm). The limitation of the matrix-type model is potentially an important issue when tracking objects with subpixel accuracy [19] or modelling long-range movements arising from frequent but small changes in an object's position [21]. The possible introduction of a three-dimensional matrix into the continuous space model would allow representation of complex cell structures (e.g. curved cell membranes, filopodia and intracellular organelles of irregular shape). Mammalian cells may contain many thousands of fluorescent molecules of one or more species [24,28]. Modelling movements, interactions and fluorescence output for a few thousand objects over a reasonable time period would require substantial computer power. It would be difficult to run such a model on a personal computer 10 years ago. However, nowadays, it possible to complete 1000-2000 image sequences in just a few minutes. For example, a computer equipped with a 3 GHz, i7 processor (Intel Co, CA, USA) runs the present model simulating 1000 cytoplasm-based molecules interacting with 1000 membrane-bound molecules and building a fluorescent image every time step at a rate of approximately 10 cycles per second. So, even at a rate of 100 fps, the model run will be only 10 times slower than real time.
The speed of the model run was increased dramatically by the algorithm calculating the probability of binding between each potential pair of reacting molecules. The solution was to calculate the AIT for the classes of molecules with known mobility and physical sizes. In most cases, the AIT rapidly decreases to negligible values at intermolecular distances greater than 0.5 mm (see figure 3 for details). Therefore, AIT values can be calculated before the model run and stored in a look-up table. This approach greatly reduces the number of calculations during the model run. There was no need to calculate which of the potential pairs was at the shortest distance-the eligible molecules would eventually bind the closest neighbour owing to the sharp increase in the AIT. This algorithm was tested by modelling PH-domain binding to the specific phospholipids on the cell membrane [14], modelling the binding of free moving membrane-bound molecules to the immobile anchors [25,27], and modelling transient dimerization of membraneembedded molecules [24]. In the case of PH-domain binding to membrane phospholipids, the dissociation rate set in the model (0.05 s 21 ) was found to be very close to the actual measured rate of 0.046 s 21 (figure 6a, and see figure 4 in [14] in order to compare it with 0.05 s 21 from the real data analysis). In the case of membrane-embedded molecules forming transient dimers (e.g. G-protein coupled receptors), the exponential distribution of dimer lifetimes produced an estimated dissociation rate, 1.2 s 21 , which was slightly higher than the rate of 1 s 21 set in the model and was very close to the measured dissociation rate (1.3 s 21 ) of muscarinic M 1 receptors studied in real two-colour imaging experiments (see figure 4 from [24] and compare it with electronic supplementary material, figure  S2). These comparisons allow us to conclude that the proposed model quantitatively reproduces the results of real single molecule experiments.
The examples presented in the Results and Discussion sections demonstrate the importance of realistic illumination and fluorescence output modelling for cross-checking single molecule experiments. It was shown that uneven illumination severely affects the conclusions of single molecule data analysis-figure 4a shows that the use of a narrow laser beam for illumination leads to the pseudo-bimodal distribution of intensities of the detected molecules and to the wrong conclusion about the presence of a mixture of monomers and dimers, whereas, in fact, the model sample contained only monomers. The summing of time steps (figure 1) is required to build realistic fluorescent images of fast-moving molecules acquired at a frame rate of 20-100 fps. Otherwise, fast-moving fluorescent molecules would look like static objects, which would affect the data analysis. This procedure is essential when different species of fluorescent molecules are tested in one experiment or when fast-moving intracellular molecules bind slow-moving or static molecules (for example, binding to cell membrane (figure 6a) or to microtubules (electronic supplementary material, figure S4f)).
This model can also be used for quantitative testing of some specific hypotheses about single molecule dynamics in cells. For example, even a modest level of permeability (10%) of a putative potential membrane barrier [2] would dramatically change the shape of the MSD-Dt plot (figure 5b) which is normally used to identify the confined diffusion [38]. Therefore, other methods should be used to probe the presence of barriers on a cell membrane. It has also been shown that single molecules moving at the membrane containing lipid rafts would have a linear shape of an averaged MSD-Dt plot (figure 5d), but the shape of the distribution of the mobility of the individual molecules would become bimodal (figure 5c, grey line) and the majority of the molecules would concentrate inside the rafts [7] if the level of mobility inside the rafts falls below 10% of the mobility outside the rafts (figure 5c, middle inset, and electronic supplementary material, video S6).
The results of modelling of some other important single molecule experiments are presented in the electronic supplementary material. These include the use of the two-colour imaging [24] for testing the monomer-dimer transitions; subunit counting in membrane-bound tetramers [36]; and modelling single molecule dynamics in the presence of intracellular structures.