OverviewThe model consists of a “simulation engine” and “model output” routines that together produce realistic fluorescence video microscopy data that is derived from known input starting values. As the model runs, complex emergent behaviors can be generated. The model system can be made as simple or complex as the experimenter requires in order to test new ideas, generate new hypotheses and validate downstream analytical approaches (i.e. image analysis methods can be tested against a ground truth). The model was written in Borland C++, but the general approach could be ported to MatLab, Python, or other computer language.Simulation engineThe approaches used for numerical modeling of diffusive and motorized movements of molecules and their interactions with each other and with structures in a simple virtual cell have been described previously26. In brief, at every time step (Δt) in the Monte Carlo simulation the following computations and logical tests are performed:
Freely diffusing molecules move a random distance in Cartesian space (x,y,z) using a Gaussian-distributed random number generator to simulate a Brownian walk. Motor-driven molecules (e.g. myosin, kinesin or dynein) move a fixed distance per time-step (Δt) along a vector to simulate motion along a cytoskeletal track (e.g. actin filament or microtubule). Many different molecules can, in principle, occupy the same voxel since their individual coordinates (x,y,z) are always held with floating-point precision.
Coordinates are checked to test if the molecule approaches within a characteristic collision distance of any other molecule and a subsequent test is made for stochastic binding probability or dissociation if the molecules are already in a bound state.
Tests are then made to check for molecular interactions within the voxel lattice elements along a linear, interpolated, path made during the time step (red circles on the arrow on Fig. 1C). The byte value of the voxel is used to lookup the voxel properties (e.g. viscosity, medium type, and properties modulating objects movements and interactions) and then suitable tests are made for stochastic binding/unbinding or reflections at voxel element boundaries. This permits interactions with cellular structures (e.g. cytoskeletal elements) and can also constrain motion of diffusing (Supplementary Material Code Example 5&6) or motorized molecules within a particular cellular compartment (e.g. cytosol, membrane or cytoskeletal filament track).
Photobleaching probability of the fluorescent tag is computed based on an exponentially-distributed random function.
The lattice structure, held in computer memory, is updated and remodeled in order to simulate cellular dynamics such as membrane protrusion, vesicle movement, cytoskeletal polymerization etc.
All of the Monte Carlo modeling computations are made and results are always stored using floating-point numbers. A rounding algorithm is then used to index the discretely sampled lattice voxels held in computer memory (memory addressing is made using a “pointer” in C + +). The byte value at the voxel address is then used to index a set of physical–chemical properties held in a 2-dimensional look-up table (i.e. a “pointer to a pointer”).Creating the 3-dimensional virtual cellComputer memory is allocated as a 3-dimensional unsigned byte array to represent “voxels” which are volume elements, analogous to “pixels” in a 2-dimensional image (Fig. 1C). The array dimensions (x,y,z) define a rectanguloid volume which encompasses the virtual cell surrounded by isotropic extracellular medium. The voxel size (scaling) is made appropriate for the purpose of the simulation; for example, a 2 nm unit length (8 nm3 voxel volume) might be used for simulating tightly curved membrane structures within a yeast cell while 10 nm unit length (1000 nm3 volume) might be sufficient for studying dynamics of molecules within the cytoplasm of a larger mammalian cell. Memory requirement in gigabytes is given by cell volume (measured in μm3) divided by the cube of voxel unit dimension (nm). So to represent a mammalian cell with dimensions (x,y,z) of 10 × 20 × 5 µm3 (= 1000 μm3) using voxels of 10 nm unit length we require 1Gbyte; whereas a higher spatial resolution of 2 nm unit length would require 125 Gigabytes of memory.Each byte of computer memory holds a code number (0–255) that is used to index a table of physical–chemical properties. For convenience values are grouped into broadly similar classes so “membrane lipids” are in the range 100–255 and “solutes” (e.g., cytoplasm, nucleoplasm, vesicular lumens, or extracellular medium) in the range 0–10 (Fig. 1C).Memory is allocated and the 3D array initialized to zero (extracellular medium) then cell compartments, cytosol, lipid membrane and cytoskeletal elements etc. are created by adding a combination of geometrical and more complex predefined shapes and using standard image processing binary morphology operations to erode or dilate in order to hollow-out or thicken the structures.Plasma membraneSimple cell shapes can be defined using a symmetrical prolate ellipsoid to represent a yeast or bacterial cell30 or a low-dimensional outline can be used to represent the more complex morphology of a mammalian cell11. The enclosed volume is initialized to represent plasma membrane (here a value of 200). A binary “erode” function is used to hollow-out the cell and fill the interior with cytoplasm (value 1) leaving a single voxel layer to represent the outer plasma membrane a subsequent binary “dilate” operation (also called a “close” operation) seals gaps in diagonal connectivity. Further dilate operations can be used to thicken the membrane if for instance the voxel unit size is < 4 nm (Fig. 1 C&D and (Supplementary Material Code Example 1&2).Sub-cellular organellesHollow spheroids, oblate and prolate ellipsoids and tubular morphologies can be added within the cytoplasm volume using Cartesian geometry and floating-point arithmetic. Bounding voxel coordinates are defined by rounding to the nearest integer value (Supplementary Material Code Example 3). To create realistic tubular networks with characteristic persistence length, Lp, ellipsoids are extruded along a variable solid angle, Δθ, at each extrusion step, Δd, using a Gaussian distributed random number generator, GRand() (with unit standard deviation) to modulate direction Eq. (1), where:$$\Delta \theta = \cos^{ – 1} \left( {e^{{\frac{ – \Delta d}{{L_{p} }}}} } \right) \times GRand()$$
(1)
Giving the average change in extrusion vector angle per extrusion step Eq. (2)$$\left\langle {\cos \Delta \theta } \right\rangle = e^{{{\raise0.7ex\hbox{${ – \Delta d}$} \!\mathord{\left/ {\vphantom {{ – \Delta d} {L_{p} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${L_{p} }$}}}}$$
(2)
In most cases, tube extrusion is terminated on contact with any non-cytoplasm voxel (see Movie S1). However, the same algorithm can be used to build filopodia projecting from the cell body into extracellular space in which case extracellular voxels are replaced with plasma membrane and cytoplasm voxels over a given extrusion distance (Fig. 1C).Voxels bounding the internal tubular and ellipsoidal structures are set to appropriate “membrane lipid” voxel values (in the range 100–255, see above) and inner volume voxels are filled with uniform “cytosol” class values (range 1–10). For further details and refinements see Code Example 4 in Supplementary Material.Seeding with fluorescently–tagged molecules: An object-orientated computing approach26,31 is used to represent each individual protein molecule of interest as an object of a certain “class” with associated properties (e.g. color of fluorescent tag, diffusion coefficient) and “methods” (or functions) that describe specific chemical or physical properties and more complex behavior e.g. motor molecules that bind and move on the cytoskeleton19,22. A random number generator was used to seed molecules at random starting coordinates within valid regions (voxels) of the cell or sub-cellular compartment (e.g. membranes, cytosol, lumen). To speed initialization, small sub-cellular compartments are filled while they are created to ensure that molecules occupy valid starting locations (see Fig. 1D). In some cases, untagged protein molecules (e.g. tag fluorescence = 0) are placed at fixed starting locations (e.g. diffusion coefficient = 0) to simulate static filamentous structures, such as actin filaments or microtubules. The principles used to build the cytoskeleton are analogous to building tubular membrane networks, described above.Dynamic cellular structuresDynamic structures (e.g. filopodia) are simulated by either adding or removing voxel elements at the tip by updating the “membrane” voxels and filling the newly enclosed or newly exposed voxel regions with cytoplasm or extracellular medium). Vesicle diffusion is simulated using a Gaussian random number generator32 to produce a Brownian walk with diffusion coefficient, D, given by Eq. (3) ref26:$$D = \frac{{k_{B} T}}{6\pi \eta r}$$
(3)
where, kBT is thermal energy, r, vesicle radius, and η, cytosolic viscosity. To reduce computation time, the (x,y,z) centroid location is held as a floating-point value and memory updates are made only when the centroid moves by > 1 voxel unit distance (see Supplementary Material, Code Example 3). To model the thermal motion of ellipsoidal organelles (e.g., mitochondria or Weibel–Palade Bodies9), rotational diffusion is computed and variation in the long-axis angle of the vesicle is updated and stored and memory voxels updated in the same way as for translational motion (above).Membrane fusion events are modeled in a probabilistic manner whenever vesicular membrane collided with plasma membrane. If the two membranes fuse, membrane proteins then diffuse between membrane compartments and molecules within the vesicle either enter or exit the extracellular space.Single molecule movements and interactionsThermal motion of individual molecules26 is computed using a Gaussian-distributed random number generator (with unit standard deviation) to give a diffusive random walk displacement along each axis (x, y, z) in Cartesian space at each simulation time step, Δt Eq. (4):$$\Delta x,\Delta y,\Delta z = \sqrt {2D\Delta t} \times GRand()$$
(4)
At each time step, the linear path taken through voxel space is checked to ensure every intervening voxel had the correct value, for example, membrane-proteins must move only through membrane voxels and path direction should follow the nearest membrane voxel path (Supplementary Material Code Examples 5&6). To ensure molecules remain correctly located within or moving on sub-cellular organelles all molecules associated with that organelle (vesicle) are moved the same discrete distance (voxel step) whenever the organelle location is updated (as described above). We use an algorithm, called “move to nearest valid voxel” (Code Example 9 in Supplementary Material), to ensure the correct localization of molecules during membrane remodeling (e.g. filopodial growth/shrinkage). In this case the location of each molecule is checked in voxel space and, if needed, the molecule is moved to the nearest voxel of the correct medium. See Movie S2 simulating diffusion of membrane molecules during vesicle fission. For simulations of single molecule diffusion within either cytosol or membrane the diffusion coefficient, D, for each molecular species is entered explicitly using experimentally or theoretically determined values.Whenever molecules collide with each other, binding and dissociation rules are applied knowing the respective physical/chemical properties26. Binding probability depends on collision distance and also a binding rate constant which are evaluated to a probability per time interval, Δt. Molecular complexes (e.g. homo-dimers, hetero-dimers or protein–ligand pairs) move as a single object, governed by the slowest diffusion coefficient and if one member is lipid bound (membrane protein) the complex then remains in membrane. Probability of complex dissociation depends on the dissociation rate constant expressed as a probability per Δt time period and tested using a uniformly-distributed random number generator (Rand())26.Multi-molecular interactions like the assembly and disassembly of protein complexes at the plasma membrane can be simulated as molecules are first allowed to aggregate and form clusters and then rapidly dissociate to monomers when the on-rate constant (binding) is set to zero. This type of simulation may be found useful to mimic spontaneous formation of signaling complexes or experimental manipulations like opto-genetic or caged-compound experiments (see Movie S3). Cytoskeletal dynamics of assembly “growth” and disassembly “catastrophe” of actin filaments and microtubules can be simulated as diffusing monomers stochastically bind/unbind at the growing/shrinking filament ends (see Movie S4).Model outputOutput data is documented either in raw binary or as a comma-delimited, ASCII file (“.csv” format) listing x,y,z, coordinates, bound state, fluorescence, and so on at every time point or as mock fluorescence video microscopy movies that closely mimic single fluorophore imaging experiments (Fig. 1B). The algorithms and equations used to generate sequences of synthetic fluorescence video microscopy movies under different illumination conditions have been described previously26. Briefly, for each video image illumination intensity is calculated for each fluorescent molecule according to its x,y,z coordinates and given illumination profile (e.g., TIRFM illumination with Gaussian beam profile). The mean number of emitted “photons” for each fluorophore is proportional to illumination intensity and photon noise is proportional to square-root of the mean value. In the simulation, each fluorophore emits photons until it stochastically photobleaches (as described above). Additionally, fluorophore “blinking” can also be simulated as stochastic fluctuations between ON/OFF states (See screenshots of model settings in Supplementary Materials). Photons emitted by individual fluorophores are centered at the appropriate x,y location on the mock image and spread over a Gaussian-shaped point-spread function (PSF) which matches the experimental imaging system (See text and Fig. S1 in Mashanov, 201426). Additional noise contributions from the camera and sample auto-fluorescence are also summed onto each mock image. The accumulated output of several model time-steps (2-to-1000) is summed onto the virtual x,y image to simulate frame integration times that are typical of video microscopy camera systems (e.g. 20 or 40 ms per frame) yielding “mock imaging rates” of 50 or 25 frames per second.Because the fluorophore locations at all time points are stored with floating point precision, output from the simulations are well-suited to DSTORM33 and other forms of super-resolution single fluorophore localization analysis. More advanced microscopy approaches like the use of polarized illumination and fluorophore dipole orientation4,34, STED, MinFlux35, fluorescence correlation spectroscopy (FCS)36 and other methods could be added to the model as different “model output” using the same underlying simulation engine. For example, FCS data can be simulated by narrowing the illumination PSF to ~ 0.2 μm (or ~ 0.05 μm for STED) and increasing the imaging rate to 10,000–100,000 fps and data can be output in ASCII format as intensity vs time, suitable for downstream auto- or cross-correlation analysis.