Article
Open access
Published: 15 May 2024

Vesicle and reaction-diffusion hybrid modeling with STEPS

Communications Biology volume 7, Article number: 573 (2024) Cite this article

1502 Accesses
3 Altmetric
Metrics details

Subjects

Abstract

Vesicles carry out many essential functions within cells through the processes of endocytosis, exocytosis, and passive and active transport. This includes transporting and delivering molecules between different parts of the cell, and storing and releasing neurotransmitters in neurons. To date, computational simulation of these key biological players has been rather limited and has not advanced at the same pace as other aspects of cell modeling, restricting the realism of computational models. We describe a general vesicle modeling tool that has been designed for wide application to a variety of cell models, implemented within our software STochastic Engine for Pathway Simulation (STEPS), a stochastic reaction-diffusion simulator that supports realistic reconstructions of cell tissue in tetrahedral meshes. The implementation is validated in an extensive test suite, parallel performance is demonstrated in a realistic synaptic bouton model, and example models are visualized in a Blender extension module.

A first-passage approach to diffusion-influenced reversible binding and its insights into nanoscale signaling at the presynapse

Article Open access 08 March 2021

NRN-EZ: an application to streamline biophysical modeling of synaptic integration using NEURON

Article Open access 10 January 2023

A minimalist model to measure interactions between proteins and synaptic vesicles

Article Open access 03 December 2020

Introduction

Vesicles play many diverse and essential roles in cell biology. For example in neurons, synaptic vesicles uptake and store neurotransmitter, releasing the contents into extracellular space upon a presynaptic signal¹, and postsynaptic endocytosis and exocytosis of vesicles containing AMPA-type glutamate receptors (AMPARs) determine expression of AMPARs in the cell surface^2,3 via the vesicular endosomal pathway⁴. Therefore vesicles play central roles both in chemical communication between neurons and in determining the strength of synaptic contacts, two vital processes of brain function. Vesicles diverse in structure, molecular composition (even within functionally similar systems across anatomies, such as in the case of synaptic endocytosis⁵), and function play many more important roles throughout cell biology and are essential to life processes.

The field of computational biology has turned its attention to modeling these key cell biology players in recent years to varying degrees of biological detail. Many of these early applications are model-specific with a focus on one particular aspect of vesicle function that is being investigated in the study often alongside an experimental observation. Examples of this include an investigation of autophagy⁶, vesicle mobility in the presynaptic space^7,8, calcium-mediated release probability on docked synaptic vesicles⁹, and the influence of cytoskeletal density on vesicle diffusivity¹⁰. A stochastic model of post-synaptic AMPAR trafficking by Tanaka and colleagues¹¹ is one of the most detailed vesicle models to date and includes many of the key features of the system such as endocytosis and AMPAR trafficking, but does have some simplifications and limitations such as drawing docking and exocytosis from random distributions instead of modeling full diffusion and with limited molecular complexity.

Whilst these early studies have paved a path to computational vesicle modeling, there is a need for a general tool or tools that have the flexibility to enable application to a variety of modeling systems. In addition, such a tool should also allow modeling of other features of the biological environment such as molecular reaction-diffusion processes and voltage-gated channels alongside the vesicle dynamics, in order to enable full modeling of the rich chemical processes behind important vesicle functions such as docking, priming and fusion¹².

In terms of designing a tool for such computational vesicle modeling, one possible approach is to model the forces that vesicles experience in the cytosolic environment explicitly by their interactions with proteins and other molecules in their environment (including water) to atomic detail, thus extending the field of molecular dynamics¹³. This approach has the benefit of high morphological and physical realism, however the drawback is that the biological time that can be achieved in simulations is severely limited on current hardware due to the high intensity of computation, and often the experimental data of such atomic detail is absent¹⁴. Although coarse-grained approaches can increase the achievable biological time somewhat to the order of nanoseconds or microseconds¹⁵, this still puts out of reach many models of interest in vesicle systems that may operate on a scale orders of magnitude higher such as synaptic vesicle endocytosis that operates on a timescale of seconds¹⁶ within a synaptic vesicle cycle operating on a timescale of minutes¹⁷.

In order to achieve these biological times, a good approach is to extend one of the current tools for subcellular modeling, bringing the many benefits of these established tools but with a vesicle-specific extension. STochastic Engine for Pathway Simulation (STEPS) is a voxel-based simulator that extends Gillespie’s Stochastic Simulation Algorithm (SSA)¹⁸ to simulate diffusion within tetrahedral meshes that are able to capture complex biological boundaries accurately¹⁹. STEPS has been extended since its initial release to support voltage calculations on the tetrahedral mesh²⁰ and parallel computation based on the Message Passing Interface (MPI) protocol²¹ via an operator-splitting framework²². The main benefit of SSA approaches is that they are fast and scalable²¹, and capture reaction-diffusion kinetics accurately as long as certain geometrical criteria are met^19,23. Although STEPS was originally designed with a focus on neural systems such as simulating molecular models of synaptic plasticity^24,25,26 or other dendritic spine molecular models²⁷, along with larger-scale dendritic systems^28,29,30, it has also been applied to other areas of biology such as to model viral RNA degradation and diffusion³¹, astrocytic calcium signaling³² and nanoparticle penetration in cancer research drug discovery³³. The power of the parallel solver within STEPS was recently demonstrated by simulation of an entire Purkinje cell at molecular detail for seconds of biological time on a supercomputer³⁴.

In this study we describe how we take the many benefits of this parallel SSA-based simulator and extend it for modeling vesicles, which cannot realistically be represented by point-molecules within the regular SSA. Our aim is to provide a flexible modeling tool that captures many important features of vesicles and their processes so that the tool can be applied to many different model systems. Vesicles are represented with enough realism to extract their key features such as mobility within cells and near membranes, their molecular interactions, and important features such as endocytosis and exocytosis, but by an approach that is capable of achieving biological times that are out of reach for molecular dynamics approaches. We describe our unique hybrid SSA-vesicle approach where vesicles and vesicle-related phenomena are given special treatment but with a firm focus on computational performance, and the benefits of the SSA are maintained for other chemical species and reaction-diffusion processes. We describe in some detail the implementation of our methods, and demonstrate biological plausibility by an extensive validation suite of all the new modeling components introduced in this paper. In addition, we discuss potential sources of error in models and how these can be avoided or minimized. We introduce the parallel MPI-based implementation, demonstrate performance in a realistic synaptic bouton model, and visualize in a Blender extension module.

Results

In this section we describe sequentially the features of vesicles that are modeled, and with results from test models by means of validating the method. We begin by qualitatively exploring vesicle mobility and crowding effects in the STEPS implementation. Later validations, however, are then generally run under well-mixed conditions (where crowding effects are negligible) in order to match to analytical solutions, although we should note in more complex models the reaction rates can be under-sampled or over-sampled due to effects of crowding.

Figure 1 shows schematically the biological phenomena that STEPS simulates in the vesicle extension. Each of these features are described within the following sections.

**Fig. 1: The full set of behaviors that STEPS implements for Vesicle and Raft modeling, described further in text.**

Vesicle overlap with tetrahedral elements in a mesh

Vesicles exist within the tetrahedral mesh environment that represents cellular volumes such as the cytosol, and vesicles may not exist outside of the mesh boundaries. Based on the position and size of the vesicles, the overlap between vesicles and the tetrahedrons is found (see Sphere-Tetrahedron overlap) and each individual vesicle as a computational unit stores information about this overlap. Each vesicle knows which tetrahedrons it overlaps, and for each tetrahedron, the volume of the tetrahedron that it overlaps, which allows all the features presented in this study to be supported such as interaction with the SSA, passive transport, and so on. As we aim to provide a general modeling tool to support a wide range of modeling scenarios there is no restriction on the size of this overlap, however Fig. 2 gives information about the range of overlap from the meshes used in the study, covering a wide range of modeling scenarios. From the largest tetrahedrons that average over 100 nm in size (Fig. 2a), synaptic vesicles of 40 nm diameter typically overlap just one or a small number of tetrahedrons, whereas Fig. 2b shows a scenario where tetrahedron size is smaller than vesicle diameter and so vesicles overlap many tetrahedrons at any given time. For the most realistic modeling scenario presented in this study (based on a realistic synaptic bouton model³⁵), tetrahedrons are of comparable size to synaptic vesicles (Fig. 2c). The modeling features are designed to work under any of these scenarios, though they can present different sources of errors, as is discussed later in Vesicle surface protein and vesicle internal molecule transport, and interaction with environment.

**Fig. 2: Size comparison between vesicles and the tetrahedral meshes used in this study.**

Vesicle mobility: passive and active transport in the crowded cell environment, partial mobility in clusters, and tethering to the cell membrane

One of the most important features of vesicles is that they exhibit various forms of mobility, dependent on many factors such as location, function and cell activity³⁶. They can, for example, be freely-diffusing within the crowded environment of the cytosol or other cell compartments (Fig. 1(i)), effectively immobile when docked to tethering molecules in the cell membrane (Fig. 1(ii)), or partially mobile when clustered to other vesicles³⁷ in functional pools³⁸ (Fig. 1(iii)). STEPS supports all of these situations. By default, vesicles freely diffuse within their cellular boundaries by Brownian motion diffusion (see Vesicle definition and diffusion). In STEPS, this is modeled as changes in position on fixed, user-defined time steps, Δt, and so the size of these diffusion steps depends on diffusion coefficient and the size of Δt. Figure 3a shows example diffusion steps for a diffusion coefficient of 0.06 μm²s⁻¹ (which is the free cytoplasmic diffusion of synaptic vesicles after hydrodynamic effects found in ref. ⁸ and also used in Fig. 3c) for typical Δt choices of 0.1 ms and 1 ms. Comparison between Fig. 3a and Fig. 2c show that vesicles usually diffuse less than the size of one tetrahedron per diffusion step under realistic modeling scenarios, although this is not a restriction on a STEPS model.

**Fig. 3: Validation of vesicle mobility, and mobility of molecules on vesicle surface.**

Because vesicles cannot overlap other vesicles, mesh boundaries, or other boundaries that might represent, for example, internal cellular structures, the mobility of vesicles is reduced in the presence of such structures⁸ thus capturing the important known features of cellular crowding³⁹. To demonstrate this, first we investigated vesicle diffusion in the case where 0%, 20%, 40% and 60% of the cell volume is occupied, showing how the software captures the expected diffusion distance in the 0% case, and mobility reduces with increasing volume occupancy (Fig. 3b), as expected. Supplementary Note 1 gives further information about this model.

In addition, we tested a mitochondria-based model introduced by Rothman et al.⁸, as detailed in Supplementary Note 2. In this model 28% of the volume is occupied by immobile mitochondria in the tetrahedral mesh, a further 25% by immobile vesicles, and 17% by mobile vesicles (the model is shown in the inset of Fig. 3c). We demonstrate that STEPS captures the expected features of this model in that the effective diffusion rate of the mobile vesicles is not constant in time; it approximates the Brownian-motion rate (0.06 μm²s⁻¹ in this model) only on short time-scales and is reduced by steric interactions on longer timescales (Fig. 3c) approaching physiological values⁸. STEPS does not model the very fast timescale effects of hydrodynamics⁸.

Vesicle active transport

Active transport on structures such as microtubules and actin filaments (Fig. 1(iv,v)) is an essential feature of vesicle mobility, especially in neurons⁴⁰. For example, synaptic vesicle trafficking along actin filaments delivers vesicles to the active zone⁴¹, and post-synaptically recycling vesicles are actively transported to the spine surface to deliver AMPARs during long-term potentiation⁴². Vesicles are transported by motor proteins myosin on actin, and kinesin and dynein on microtubules⁴³. A feature of the motor-protein walk is static periods punctuated by sharp steps of regular size to the next position along the structure^44,45,46. The dwelltimes of which are near single-exponential or double-exponential^45,47,48, indicating one or sometimes two rate-limiting steps.

STEPS models active transport by a system of branched Paths representing cytoskeletal filaments that vesicles may interact with if they cross them in their environment. Vesicles travel along the path from the beginning point to the end point by a user-defined speed and stepsize (for example the stepsize of myosin V is 36 nm, consisting of a 25 nm working stroke and an 11 nm diffusive movement⁴⁹). The default behavior is to sample dwelltimes from a single-exponential distribution, although a double-exponential may optionally be used. The total distance traveled divided by the sum of the dwelltimes then equals the given speed, with noise arising from the distribution of dwelltimes (instead of a single value). Supplementary Note 3 and Supplementary Fig. 1 demonstrate the implementation of single-exponential and double-exponential dwelltimes in STEPS and Supplementary Video 1 shows vesicles undergoing active transport on a path in STEPS in a simple model.

As well as active transport, paths may be useful in other modeling contexts such as vesicle binding to actin cytoskeleton in cluster formation⁵⁰ and active zone tethering⁵¹.

We validate the implementation of paths in a model measuring vesicle displacement on path vs expected displacement, as demonstrated in Fig. 3d. The model vesicles undergo several periods of active transport along paths at the expected speeds (in this example with small stepsize to counteract stochastic effects for the purpose of validation). The implementation of this model is described in further detail in Supplementary Note 4.

Vesicle surface molecule diffusion

Vesicles contain proteins in their surface which may have been present at the site of endocytosis (Fig. 1(vi)) and which can undergo various interactions such as phosphorylation upon interaction with the cytosolic environment. These vesicle surface molecules may play important roles such as neurotransmitter uptake (Fig. 1(vii)), docking to cell membrane (Fig. 1(viii)) or binding to other vesicles (Fig. 1(iii)). These molecules may be mobile in the vesicle surface, undergoing lateral diffusion (Fig. 1(ix)). This amounts to modeling diffusion on a spherical surface, by which STEPS applies the efficient approximate propagator algorithm of Ghosh et al.⁵². As described further in Supplementary Note 5 and in Supplementary Fig. 2, to test the STEPS implementation of the Ghosh algorithm we tested a regular diffusion model on a spherical surface mesh in STEPS and compared to the Ghosh algorithm (Fig. 3e shows one example of angular displacement after time 2 ms of a surface molecule diffusing at 1 μm²s⁻¹). The Ghosh algorithm compares accurately to the measured angular displacement distribution from our surface diffusion model (Fig. 3f) and is a reliable, fast algorithm for implementing diffusion of surface molecules on vesicles in STEPS.

Exocytosis, endocytosis, fusion and budding

Exocytosis (Fig. 1(x)) is an essential cell process by which luminal species inside vesicles can be secreted into the extracellular space (for example neurotransmitter release at chemical synapses) and membrane proteins may be inserted into the cell membrane. Both processes are vital to cell function. Vesicles are also able to fuse to internal membranes such as endosomes (Fig. 1(xi)), and both this and exocytosis are modeled in STEPS by the same computational framework within the SSA. STEPS is able to model the complex chemical processes leading up to and resulting in exocytosis, such as those involved in docking, priming and fusion (Fig. 1(viii))⁵³.

For synaptic vesicles it has been proposed that kiss-and-run fusion⁵⁴, whereby a vesicle does not completely collapse upon neurotransmitter release but remains intact and is immediately endocytosed and recycled, is utilized in neurons due to the clear advantage of speed of recycling. Full-collapse fusion and kiss-and-run are proposed to coexist in many systems, although the prevalence and significance of kiss-and-run fusion is somewhat debated (for reviews see refs. ^55,56). In addition, another mode of release known as open and closed fusion has been described and, like kiss-and-run, does not result in full collapse of the vesicle but is distinct from kiss-and-run for example by amount of vesicle content released and is proposed as the main mode of exocytosis⁵⁷. Partial release of luminal content during open and closed fusion has been directly measured, for example Li et al.⁵⁸ measure 64% of total catecholamine content release during stimulated exocytosis in PC12 cells. Since a modeler may want to employ any of these modes of exocytosis, STEPS supports modeling of full-collapse fusion, fusion where the vesicle is preserved (that may be used to model kiss-and-run or open and closed fusion), and partial release, with a user able to employ any of these different modes within the same model.

Within a full-collapse fusion exocytosis event in STEPS all luminal molecules are released into the extracellular space (Fig. 1(xii)) and all vesicle surface-bound species are transferred to the cell membrane (Fig. 1(xiii)). When fusion is not full-collapse and the vesicle is preserved vesicles can partially preserve luminal composition⁵⁵, and so in STEPS all or partial luminal contents are released (Fig. 1(xii)) at the modeler’s discretion. As mentioned, fusion is also possible on internal structures (Fig. 1(xi)), enabling modeling of cycles such as the the post-synaptic AMPAR cycle by which AMPAR-containing vesicles fuse to endosome⁴.

Figure 4 a provides a basic validation of exocytosis in STEPS within the SSA framework, with further details of this model given in Supplementary Note 6. Figure 4e visualizes a full-collapse fusion event in STEPS of one synaptic vesicle in the active zone, where release of neurotransmitter into the extracellular space and retention of the SNARE complex and vesicle surface-bound species in the cell membrane can be seen. Supplementary Video 2 shows multiple exocytosis events in the same active zone model.

**Fig. 4: Validation of Endocytosis and Exocytosis in STEPS by basic test models where analytical comparisons can be made.**

Endocytosis (Fig. 1(vi)) from an invaginated region of a cell membrane, and budding from internal structures such as liposomes⁵⁹ (Fig. 1(xi)) both involve the creation of vesicles, and are modeled in STEPS as events within the SSA. These events can optionally be modeled with a species dependency, meaning that the event is only active and available for application within the SSA if a certain species signature is met, for example if a certain number of clathrin molecules have been recruited to the region⁶⁰ defined as an Endocytic Zone within STEPS (Fig. 1(xiv)). This region may be defined as a collection of surface triangles in the tetrahedral mesh, or as molecular population of a Raft (Rafts are described later in section Membrane Rafts: Formation, diffusion and interaction with their environment). Figure 4b–d provides a basic validation of endocytosis in STEPS within the SSA framework when the endocytosis region is either based on rafts (Fig. 4b) or comprised of mesh surface triangles (Fig. 4c).

STEPS is capable of modeling the complex molecular machinery that generates endocytosis events, but does not yet support adaptive meshes and so the geometric changes are not fully modeled. Instead, when an endocytosis event is chosen to take place, the spherical space just inside the endocytosis region is reserved and a vesicle is formed in place. This vesicle will contain the molecular makeup of the zone in which endocytosis occurred in its membrane, to mimic the molecular biology (Fig. 1(vi)). Figure 4d validates the expected production of vesicles during repeated endocytosis events. As dependent species are consumed, the rate of production of vesicles can be seen to decrease in this validation model. The endocytosis validation models are described in more detail in Supplementary Note 7.

Vesicle surface protein and vesicle internal molecule transport, and interaction with environment

Vesicles contain proteins embedded in their membrane which allow them to interact with their environment through protein-protein or other interactions (Fig. 1(vii, viii)), or even interact with each other (Fig. 1(iii)). An important example of such interactions are those that allow synaptic vesicles to go through the process of docking and priming, which is controlled by interactions between vesicle proteins and proteins embedded in the cell membrane⁶¹ (Fig. 1(viii)). Further, primed vesicles then sense calcium released into the cytosol to undergo the process of fusion and neurotransmitter release through a cascade of chemical interactions^53,62.

STEPS supports all of these situations by addition of a new type of reaction termed a Vesicle Surface Reaction by which reactants may exist in a number of regions: they may be compartment-based (e.g. cytosolic) molecules, molecules in the vesicle surface, or molecules embedded in the cell membrane or other membranes. A reaction between a molecule in the vesicle surface and the cell membrane can be used to model docking (Fig. 1(viii)), for example. In addition, reaction products may appear internally within the vesicle to model e.g. neurotransmitter pumps and filling (Fig. 1(vii)). These reactions are then solved stochastically as regular SSA reactions. This allows modeling of the rich interactions that vesicles may undergo with their environment and allows, for example, the entire synaptic vesicle cycle to be modeled to full molecular detail³⁵. To faithfully model docking reactions specifically, optionally a maximum distance may be defined so that the reaction will only occur within a specified distance of the membrane surface, and vesicles may be immobilized when effectively docked close to the membrane (Fig. 1(ii)). The reverse reaction can also be modeled where vesicles may become undocked and mobilized, which is important if modeling kiss-and-run or open and closed fusion.

Figure 5 demonstrates our validation of the vesicle surface reactions by way of interaction with vesicle surface molecules and cytosolic species. We extend our regular validations¹⁹ for these interactions, and demonstrate very close agreement in four different models, namely first order irreversible (Fig. 5a) and reversible (Fig. 5b) reactions, and second-order irreversible (Fig. 5c) and reversible (Fig. 5d) reactions, in every model demonstrating very close agreement to the expected solution analytically. Supplementary Note 8 contains further details of the model setup and parameters used.

**Fig. 5: Validation of vesicle surface molecule interactions.**

Vesicle interactions with each other

A known behavior of vesicles is their ability to form functional clusters or pools^38,63 in which mobility is strongly reduced⁶⁴. It has been proposed that such clusters are formed by vesicles tethering to each other for example in the case of synaptic vesicles by synapsin interactions^50,65. STEPS models this phenomenon by allowing vesicles to loosely bind to each other through interactions between their surface proteins. A special type of reaction termed a Vesicle Binding event models this, alongside the reverse process Vesicle Unbinding (Fig. 1(iii)).

As vesicles bind in STEPS they form a computationally distinct complex termed a Link Species on their surface (Fig. 5f). Link Species have the property of length, which is by a user-defined upper and lower boundary, and the Vesicle Binding event can only occur if and when the resulting Link Species formed would be within its permissible bounds. Link Species allow some mobility of linked vesicles as long as any movement results in the Link Species remaining within its bounds. In addition, Link Species can undergo all the usual chemical interactions such as interaction with cytosolic or vesicle-bound molecules, and diffusion on the vesicle surface. Through this behavior, vesicles are therefore able to form clusters in which they exhibit loose, partial mobility, phenomenologically similar to their low mobility in pools⁶⁴. Thus vesicle clustering and pool formation such as the classical reserve pool, recycling pool and readily releasable pool³⁸ can be faithfully modeled in STEPS, along with their interactions such as activation by calcium, and the process can be modeled to high molecular detail.

We provide a validation model of the Vesicle Binding reaction in STEPS, described in Supplementary Note 9, which has a known analytical solution to which STEPS output can be compared. Figure 5e shows close agreement between the STEPS model and analytically expected behavior of this model. In addition, Fig. 5f demonstrates cluster formation in which synapsin dimers bind vesicles to each other, showing liquid phase behavior³⁷.

In addition, Supplementary Video 3 shows a cluster model in which vesicles are bound to each other by synapsin dimers yet still exhibit some mobility. The video shows 300 vesicles that are bound together and 10 inert vesicles that do not interact with the cluster but can be seen to mix with it, demonstrating liquid phase properties of the cluster.

Membrane rafts: formation, diffusion and interaction with their environment

Nanoscale domains of clustered proteins, lipids and other molecules in cell membranes play many important roles in signaling and trafficking⁶⁶. Referred to as Lipid Rafts, these important membrane subdomains span a large range of sizes, reported as 10 nm to 200 nm in a definition that includes caveola⁶⁷, although a recent fluorescence resonance energy transfer study measured a lower bound of 5 nm⁶⁸. Lipid rafts are dynamic in their molecular composition, and may be mobile and diffuse freely in the membrane or become bound to cytoskeleton and immobile⁶⁹. Somewhat analogously to vesicles within cytosol and other cell compartments, they cannot be modeled as regular volumeless molecules within SSA methods and require special treatment.

We implement Rafts in STEPS as regions that occupy an exclusive, fixed radius within cell membranes. To mimic the biology, rafts may be mobile by a user-defined diffusion rate within cell surfaces (Fig. 1(xiii)) or immobile. They can undergo all the usual chemical interactions with their environment including with cytoplasmic molecules and surface-bound molecules and, in addition, rafts can be created upon exocytosis of a vesicle (as proposed for example in ref. ⁵¹) or form an endocytic zone from which a vesicle will be produced, sharing the surface molecular composition with the exocytosed or endocytosed vesicle (Fig. 1(vi, xiii)).

Since rafts are associated with a triangular meshed surface in a STEPS simulation, we apply the same diffusion rule for rafts as for regular SSA surface species. Figure 6a validates this implementation in a diffusion model where crowding effects are minimal, although it should be noted that diffusion rates can be affected by crowding if raft density is high. In addition, and similarly to our validation of vesicle surface reactions, we validate our implementation of raft surface reactions by models of first-order irreversible (Fig. 6b) and reversible (Fig. 6c) reactions, second-order irreversible (Fig. 6e) and reversible (Fig. 6f) reactions. As a useful modeling feature, it is also possible in STEPS to model formation of rafts when a particular species signature is met, that is when a user-defined list of species are present at a given number within a triangle. Conversely, dissociation and dispersion can be modeled if a raft molecule composition goes below a certain threshold. Figure 6d validates these implementations in a simple model that includes both of these phenomena. Further details of all of these models can be found in Supplementary Notes 10 and 11.

**Fig. 6: Validation of rafts and their interactions.**

Computational performance in parallel

STEPS implements a parallel splitting algorithm for the vesicle solver, since an original serial implementation was found to be very limited in applications due to all the extra computation involved in supporting vesicle and raft modeling. Described in further detail in MPI implementation, we achieved an approximately 50X speedup on a realistic and detailed synaptic bouton model of Gallimore et al.³⁵, simulating 1 ms of biological time in 7.3 s of wall-clock time (averaged over a total run of 2 s biological time). This performance is sufficient to bring biological systems that operate on timescales of seconds or minutes such as the synaptic vesicle cycle into scope. If computer cluster runtime constraints are a concern, checkpointing further extends the biological time that is achievable. Further improvements to MPI performance may be possible by parallelizing the vesicle routines and with further optimization to the MPI communication schemes in the future.

Discussion

In this study we presented a new modeling tool that supports realistic simulation of vesicles and their many important biological functions such as uptake, transport and secretion, as well as lipid rafts. To our knowledge, this is the first approach to successfully combine such modeling in a hybrid framework with spatial SSA processes, enabling detailed chemical modeling of important biological systems such as the synaptic vesicle cycle. Every stage of development was carefully tested and validated, with many of the validations described in this paper, and implemented in realistic biological scenarios. In one example described by Gallimore at al.³⁵ and shown in MPI implementation, a detailed pre-synaptic bouton model simulated the synaptic vesicle cycle to full molecular detail.

The software was implemented in an MPI framework that took into account the unique challenges of modeling volumetric objects within the partitioned spatial SSA framework. Although this is a difficult system to parallelize, by taking advantage of the fact that vesicles diffuse much slower than smaller molecules we were able to sufficiently improve performance by orders of magnitude compared to a serial implementation, putting into reach biological systems operating on timescales of seconds or minutes. Such timescales are unreachable in molecular dynamic approaches.

We have developed this as a general tool that has the flexibility to enable application to a variety of modeling systems with realistic modeling of vesicles and lipid rafts and with the power to reach size and time scales of interest. The software is free for anyone to use under the GNU General Public License v3.

Methods

The new entities and kinetic processes in STEPS 5.0, and the relationship between them

STEPS follows an object-oriented programming style, and many new computational objects have been added in 5.0 for vesicle modeling, including Vesicles, Rafts, and their modes of interaction with their environment such as Vesicle Surface Reactions, Exocytosis and so on, as described in Results. Figure 7 summarizes the relationship between these objects with a color-coded description, and Supplementary Table 1 shows this information in tabulated form. One important feature is that molecules on a vesicle surface (both Species and Link Species) have an exact position so that they can both undergo diffusion on the surface of the vesicle (Vesicle surface molecule diffusion), and so that their position can be found within the mesh for availability to the SSA (for example for a Vesicle Surface reaction or Vesicle Binding event). Since there is no concept of diffusion on a Raft surface in STEPS, molecules contained in Rafts do not have an exact position but are randomly assigned to Raft overlap triangles for availability to SSA processes (such as Raft Surface reactions).

**Fig. 7: Schematic summary of kinetic processes.**

Vesicle definition and diffusion

Vesicles are approximated in STEPS as spherical volumes as obviously the most natural choice of the available geometrical primitives. This volumetric modeling of vesicles is a clear difference to the effectively point-particle SSA species in STEPS used to model smaller molecules such as ions, ion-channels and enzymes. Cellular vesicles vary in size; for example in neurons synaptic vesicles are approximately 40 nm in diameter and can transport a relatively small amount of neurotransmitter and surface proteins but can release quickly, whereas dense-core vesicles that transport neuropeptides are larger at 80 nm to 200 nm and operate more slowly⁷⁰. A STEPS user specifies vesicle diameter that vesicles of this type will then assume during the simulation as effectively a sphere that operates within the realms of the tetrahedral mesh environment. A STEPS user is not limited in how many different types of vesicles they may define and so, for example, clear synaptic vesicles and dense core vesicles can both be modeled within the same simulation with differing properties, such as size, rate of diffusion, fusion time course, membrane protein composition, and so on. The user may then also define a fixed timestep on which to update vesicle positions by diffusion, although this will default to 1 ms if not specified.

In order to model Brownian motion on these regular timesteps, Δt, a diffusing vesicle’s next position is chosen by sampling the x, y, z displacement separately:

$$\begin{array}{rcl}\delta x&=&r1* \sqrt{2D{{\Delta }}t}\\ \delta y&=&r2* \sqrt{2D{{\Delta }}t}\\ \delta z&=&r3* \sqrt{2D{{\Delta }}t}\end{array}$$

where r1, r2, r3 are random numbers generated on the standard normal distribution. The new position is then tested to see if it is a valid position or not. For example, the vesicle must not overlap any other vesicle nor must it cross a boundary such as the mesh surface. To optimize the search of other vesicles, only vesicles known to overlap the tested tetrahedrons are checked. To calculate whether a vesicle crossed a compartmental boundary requires knowledge of the sphere-tetrahedron overlap of the new position (see Sphere-Tetrahedron overlap).

Figure 3 a shows diffusion distances for a diffusion coefficient of mobile synaptic vesicles of 0.06 μm²s⁻¹ ⁸ for timesteps of 0.1 ms and 1 ms. Since diffusion distance is lower than vesicle diameter, meaning that the vesicle will continuously overlap some volume when diffusing, these appear to be reasonable choices for this type of vesicle, allowing a more optimal overlap search (see Sphere-Tetrahedron overlap). Sources of error, and reporting to user expands on when the timestep can be too large, and how this is reported to the user, with suggestions for a better value.

Sphere-Tetrahedron overlap

STEPS uses an external library to determine the sphere-tetrahedron overlap⁷¹ so that the set of tetrahedrons that a vesicle overlaps can be found. For application of this algorithm, it is essential to avoid testing all tetrahedrons in the entire mesh due to runtime concerns. If a vesicle is diffusing by a diffusion distance that is smaller than the vesicle’s diameter, the set of tetrahedrons previously overlapped are taken as the starting point for the overlap search because there is guaranteed overlap in this set. If the diffusion distance is bigger than the vesicle diameter, a walking algorithm is used to find a path between the current vesicle center and the future vesicle position. The algorithm consists in a variant of the A* algorithm⁷²: a set of tetrahedrons is initialized with the tetrahedron containing the current center of the vesicle, neighboring tetrahedrons are then explored by processing first tetrahedrons closest to the destination. The algorithm stops when the tetrahedron containing the destination has been found or when all tetrahedrons that are below a threshold distance ${d}_{\theta }=3\parallel \overrightarrow{{p}_{s}}-\overrightarrow{{p}_{e}}\parallel$ from the destination have been explored (with $\overrightarrow{{p}_{s}}$ the starting position and $\overrightarrow{{p}_{e}}$ the destination). With this distance threshold, this algorithm prevents vesicles from jumping over gaps in the mesh. When the vesicle diameter is small compared to the size of a tetrahedron, it is possible for d_θ to be too small and prevent exploration of even the immediate tetrahedron neighbors. To avoid this, the algorithm always explores at least 10 tetrahedrons before stopping.

A breadth-first search is then used to test overlap of tetrahedrons, stopping if no overlap can be found within a layer or overlap reaches 100% of the vesicle volume (which means that the new position is good). Overlap under 100% means that the vesicle is not full covered by tetrahedrons and so the position being tested is outside of the permitted environment (i.e outside of compartmental boundaries).

Vesicle-mesh overlap effects on molecular reaction and diffusion

To approximate the effect of vesicle presence on the mobility of other molecules (modeled within the operator-splitting framework²²), we apply a correction to the finite-volume diffusion rate¹⁹ of molecules proportional to the volume occupancy. The modified diffusion rate from tetrahedron k to tetrahedron l of species i is:

$${d}_{k,l}^{* }=\frac{{D}_{i}{A}_{k,l}}{{V}_{k}^{* }\delta {x}_{k,l}}\,$$

where D_i is the macroscopic diffusion coefficient of species i, A_k,l is the face area between tetrahedrons k and l, ${V}_{k}^{* }$ is the reduced volume by vesicle occupancy, and δx_k,l is the barycenter-barycenter distance. This reduced volume effect intuitively reduces the dwelltime for a molecule inside a tetrahedron proportional to the reduced volume, whilst also ensuring SSA species cannot occupy a fully overlapped tetrahedron and so do not occupy the same space as a vesicle. Note that mathematically the diffusion rate will approach infinity as the reduced tetrahedral volume approaches zero (meaning complete coverage by a vesicle) and is effectively replaced by a high rate for diffusion to a neighbor within the operator-splitting framework if this situation occurs. In realistic meshes, however, usually full coverage does not occur (Fig. 2c) (and won’t occur unless tetrahedrons are very small, Fig. 2b) and in addition molecules are expected to evacuate tetrahedrons before full overlap occurs. Diffusion to a fully-overlapped tetrahedron is not permitted.

We tested the effect on regular diffusing SSA species in a model with 10% volume occupancy by vesicles, as described in Supplementary Note 12. In all cases, the tetrahedrons captured the binomial distribution of molecules with low error, as shown in Supplementary Fig. 4. The error in the fit to the binomial distribution was always less than 0.5%. This means that no spatial bias is introduced on the SSA species distributions by the modified diffusion rates.

STEPS also applies a volume correction to SSA reaction rates in order to capture effects of a reduced volume. For example, the rate, c, of a particular second-order reaction taking place in tetrahedron j is:

$${c}^{* }=\frac{K}{{N}_{A}{V}_{j}^{* }}\,$$

where K is the macroscopic reaction rate, N_A is Avogadro’s number and ${V}_{k}^{* }$ is the reduced volume by vesicle overlap. In this sense, the volume that vesicles occupy is not regarded as part of the well-mixed volume of a tetrahedron. This can be expected to accentuate small-volume errors on reaction rates⁷³ if overlap is high. However, in practice, because overlap is transient and often relatively small, this effect also may be small or negligible because the effective volume is still usually within the accurate size window¹⁹. It should be noted that all validations presented in this paper, such as the vesicle surface molecule interactions presented in Vesicle surface protein and vesicle internal molecule transport, and interaction with environment and Vesicle interactions with each other pass within this framework. Unimolecular reactions are unaffected. The following section describes under which circumstances the second order reaction error can become large, how this is reported to the user, and steps that can be taken to reduce the error.

Sources of error, and reporting to user

In this section we describe sources of error that can appear in models, often either by model design or from mesh construction, and how STEPS reports these expected errors and suggests modifications to the model in order to reduce or nullify these errors.

Basic checks on STEPS models

Firstly, when a STEPS simulation is created, a series of automatic checks are run on the model parameters to flag any potential mistakes or problems, with warning messages displayed to the user. For example, if any molecular species only appear on the right-hand side of reactions then a warning is displayed because nothing in the model consumes the species, and reaction rates are compared to other reactions in the model of similar type to check for any obvious outliers that could be simple typing errors. We have extended these checks for the new reactions that take place in STEPS 5.0. In addition, basic checks specific to vesicle and raft modeling have been added. STEPS checks whether the sizes of vesicles and rafts are larger than the size of the mesh and whether diffusion coefficient for vesicle surface diffusion are so high that one vesicle Δt would result in a random position on the surface, since these are indicative of modeling errors.

We now describe more detailed checks on diffusion and reaction rates, and when warnings and suggestions are displayed to the user.

Vesicle diffusion

Perhaps the most obvious source of error when modeling vesicle diffusion on fixed time steps, Δt, is if the usual diffusion distance is too large. We suggest that this is the case when a high proportion of diffusion steps result in a vesicle diffusing greater than its diameter plus the smallest diameter of other vesicles in the model. A diffusion step of this distance will, in theory, allow a vesicle to pass through other vesicles whereas they would meet these obstacles if diffusing with a lower distance, and therefore a large diffusion distance could underestimate effects of crowding. Figure 8a shows, in the case where all vesicles in the system are of the same diameter, the proportion of diffusion steps over 2 diameters, as a function of the root mean squared diffusion distance. The fraction of vesicle diffusion steps above a distance d can be obtained from $1-F(\frac{3{d}^{2}}{6D{{\Delta }}t};3)$ with F(x; k) the cumulative distribution function of the chi-square distribution with k degrees of freedom and D the vesicle diffusion coefficient. STEPS checks the model parameters and reports to the user if 5% of the diffusion steps are expected to be over a vesicle’s own diameter plus the smallest diameter of other vesicles in the system, and suggests a vesicle Δt to the user at which this ceases to be the case.

**Fig. 8: Sources of error in vesicle models in STEPS.**

Bimolecular reactions

As previously noted, and described elsewhere^19,73, tetrahedron size can influence the simulated reaction rate of bimolecular reactions where molecules must diffuse to meet each other before reacting. Although, for realistic parameters, a voxel size window usually exists where relative error is low^19,23, the error can become large if tetrahedron sizes go below a certain size threshold- an important consideration for mesh construction.

The new components in STEPS 5.0 that support bimolecular reactions modeling and so are subject to these tetrahedron size considerations are:

Vesicle Surface reaction (bimolecular)
Raft Surface reaction (bimolecular)
Vesicle Binding

We investigate this effect with a Vesicle Surface Reaction model where reactions take place between reactant molecules in the cytosol and on the vesicle surface. In Fig. 8b we demonstrate the error in our simulated reaction rate of the vesicle surface reaction validation model (Vesicle surface protein and vesicle internal molecule transport and interaction with environment), but with different meshes representing the same morphology of a 0.5 μm diameter sphere. 8 meshes were used with average tetrahedron size ranging from 120 nm (291 tetrahedrons) to 12.8 nm (265,307 tetrahedrons) (these two particular meshes can be seen in Fig. 2a, b). This spans approximately an order of magnitude in size, and therefore 3 orders of magnitude in volume. We simulate two cases, one where the vesicles are diffusing at 0.1 μm²s⁻¹ and one where vesicles are immobile, which could represent docked or clustered vesicles. We observe that the error generally stays low until 12 nm, although the error is usually slightly higher with immobile vesicles until this point. At 12 nm the error in the reaction rate with immobile vesicles increases sharply to 6.05% whereas it remains low with the diffusing vesicles, at 2.78%. This is because with immobile vesicles, the relative diffusion rate between the reactants is reduced, amplifying the small volume effect. Vesicle modeling can contain immobile or low-mobility vesicles by effects such as docking or clustering, and we conclude that it may be more difficult to capture bimolecular reaction rates under these scenarios with very small tetrahedrons. For these reasons, and for the observations in Fig. 8b and our previous observations¹⁹, as well as observing the typical sizes of tetrahedrons in a realistic biological mesh as shown in Fig. 2c, STEPS checks a mesh and reports to the user if 10 percent of tetrahedrons are below a size of 20 nm. This check also captures simple mesh scaling errors.

Since we observed that the vesicle diffusion rate can influence bimolecular reaction rates, we can expect that the vesicle diffusion timestep may also have an effect. We investigated this under the scenario of immobile cytosolic reactants as a tractable mathematical problem and effectively worst case modeling scenario, where the observed error is smaller with mobile cytosolic reactants (see Supplementary Fig. 5 for cytosolic reactants diffusing at 1 μm²s⁻¹). Under this scenario, bimolecular reaction rates can be underestimated if the vesicle Δt is too high. For higher vesicle Δt, vesicle surface species spend a longer time in each tetrahedron than would be expected if vesicle diffusion was continuous, immobile reactants get more depleted than they should and this can lead to an underestimation of bimolecular reaction rates. Figure 8c shows the error on the estimated rate as a function of vesicle Δt and tetrahedron sizes. Based on the sizes of tetrahedrons in the mesh, STEPS computes a threshold value for vesicle Δt (red line) above which the error on reaction rate is deemed too high, and displays a warning to the user, suggesting a more appropriate vesicle Δt. This threshold is computed with ${{\Delta }}{t}_{\theta }=\frac{{r}_{tet}^{2}}{15D}$ which corresponds to the average time needed for a point diffusing in 3D with diffusion coefficient D to exit a sphere of radius r_tet when its initial position is uniformly distributed in the sphere. For any given mesh r_tet is taken so that the volume of the sphere is equal to the average volume of a tetrahedron. Details related to how the threshold was derived are available in Supplementary Note 13 and Supplementary Eqs. (10)–(19).

We should note that some reactions are not affected by the tetrahedron size consideration because they do not involve two reactant molecules. In STEPS 5.0 these are:

Vesicle Surface reaction (unimolecular)
Raft Surface reaction (unimolecular)
Vesicle Unbinding
Exocytosis
Endocytosis
Raft endocytosis
Raft generation
Raft dissolution

As an example, Supplementary Fig. 3 shows how a unimolecular vesicle surface reaction is not affected by tetrahedron size.

MPI implementation

We previously parallelized the reaction-diffusion and voltage solution in STEPS²¹, yet the complex additional vesicle-related functionality requires a customized solution that builds on this scheme. Our approach is to implement a splitting method whereby rank 0 (which we term the VesRaft rank) carries out vesicle-related computations such as vesicle and raft diffusion, diffusion of vesicle surface species, and application of endocytosis and exocytosis (Fig. 9a, b), and the other available processes (which we term the RDEF ranks) carry out the regular SSA calculations based on our previous method, which now however include many new phenomena such as reactions on the vesicle and raft surfaces, fusion and so on (Fig. 9a, b). This is achieved by effectively partitioning vesicles and rafts on the RDEF cores into objects that we term vesicle and raft proxies (Fig. 9a). This is essential to maintain the partitioned-SSA method that is known to scale well in realistic reaction-diffusion models²¹, with many vesicle-related phenomena operating within the SSA in parallel. Of course this requires frequent communication between the VesRaft core and the RDEF cores (Fig. 9b).

The new solver inherits the MPI communication scheme for reaction-diffusion and the voltage solutions in the previous operator-splitting implementation²¹, while adding communications of complex data structures between the VesRaft rank and the RDEF ranks (Fig. 9b). These two types of solver contain different local data, thus extra synchronization in data inquiry functions are also needed. These new challenges are addressed by the implementation of customized MPI datatype encapsulation and conditional MPI operation templates. To simplify and reduce the usage of MPI communication of multiple data entries in simulation objects (vesicles for example), data entries in object classes are firstly defined and encapsulated as custom datatypes using MPI’s user-defined datatype mechanism. Such encapsulation, together with our STL (Standard Template Library) container based MPI broadcasting and gathering templates, allows us to pack objects with multiple primitive data types as a single vector and synchronize using minimal MPI calls. This helps to reduce the latency due to frequent synchronization of small size data entries and improves the simulation efficiency. Data output in the VesRaft rank is often irrelevant to RDEF ranks, and vice versa, therefore global synchronization is not necessary, particularly for vesicle-related data. We allow the user to switch communication mode for data inquiry functions by setting two condition flags in the solver: first, whether global output synchronization should be enabled; and second, if global output synchronization is disabled, which rank should receive the actual output. Note that network communication may still occur if the acquired data is not available in the output rank solver. In practice, we enable global synchronization by default, but disable it and set the VesRaft rank as output rank for heavy vesicle manipulation routines. Internally, this functionality is accomplished by a series of conditional MPI communication and data reduction templates, which automatically identify the data availability and perform the necessary data synchronization according to the conditions.

Figure 9c shows performance of our MPI-based solution on the model of Gallimore et al.³⁵ (visualized in Fig. 9d) run for 2 s biological time on a computer cluster with varying numbers of cores. Simulation wall-clock time reduces from 386 s on 2 cores (the lowest number of cores possible in our implementation) to 7.3 s on 512 cores per 1 ms of biological time, an approximately 50X speedup.

Visualization with Blender

With the release of STEPS 4.1, the possibility to save simulation data to the HDF5 format was added. This data can then be loaded into scientific data visualization software such as Paraview⁷⁴. Although well suited to the visualization of purely mesh-based data, this approach was not ideal for the hybrid SSA-vesicle simulations described in this paper. We therefore developed a separate stepsblender python package that can read simulation and mesh data automatically saved to an HDF5 file and visualize it in the Blender 3D computer graphics software (http://www.blender.org). Simulations can be visualized interactively or rendered as movies directly from Blender. The visualization includes vesicles, rafts, link species, regular species, as well as compartments, patches, vesicle paths, and endocytic zones. Fig. 9d shows an example still shot from the model described by Gallimore et al.³⁵ and output from the Blender extension module is also shown in Figs. 3b, 4e and 5f.

This new python package contains a command line tool to import the data into Blender as well as python modules that can be imported in python scripts to further customize the appearance of model objects. To avoid memory issues in Blender, only the data corresponding to the currently visualized timestep is loaded into Blender. The data loading from HDF5 file is handled in a separate process, potentially running on a different machine, allowing the visualization of data from a computing cluster without downloading the full HDF5 files.

Statistics and reproducibility

Detailed information of the simulations shown in this study is provided in Supplementary Notes 1–13, including how many times each simulation was repeated with different sequences of random numbers, for which the mean results of these simulations are shown in Figs. 3–6 and 8. Where standard deviation is shown, this is calculated with Python package NumPy from the sample size indicated in figure legends. All results shown in this manuscript, apart from the parallel performance study (shown in Fig. 9c), may be reproduced by model scripts available at https://github.com/CNS-OIST/STEPS_Validation/tree/main/vesicles, which also includes instructions for how to run the simulations, collect data and plot the results. The script for the parallel performance study is available on request.

Data availability

All data used in this study are available in Zenodo at: https://doi.org/10.5281/zenodo.10901994⁷⁵.

Code availability

STEPS 5.0.1 is available from the STEPS public release Github repository at: https://github.com/CNS-OIST/STEPS/releases/tag/5.0.1⁷⁶.

References

Chanaday, N. L., Cousin, M. A., Milosevic, I., Watanabe, S. & Morgan, J. R. The synaptic vesicle cycle revisited: new insights into the modes and mechanisms. J. Neurosci. 39, 8209–8216 (2019).
Article CAS PubMed PubMed Central Google Scholar
Lledo, P.-M., Zhang, X., Südhof, T. C., Malenka, R. C. & Nicoll, R. A. Postsynaptic membrane fusion and long-term potentiation. Science 279, 399–403 (1998).
Article CAS PubMed Google Scholar
Lüscher, C. et al. Role of ampa receptor cycling in synaptic transmission and plasticity. Neuron 24, 649–658 (1999).
Article PubMed Google Scholar
Parkinson, G. T. & Hanley, J. G. Mechanisms of ampa receptor endosomal sorting. Front. Mol. Neurosci. 11, 440 (2018).
Article CAS PubMed PubMed Central Google Scholar
Gan, Q. & Watanabe, S. Synaptic vesicle endocytosis in different model systems. Front. Cell. Neurosci. 12, 171 (2018).
Article PubMed PubMed Central Google Scholar
Martin, K. R. et al. Computational model for autophagic vesicle dynamics in single cells. Autophagy 9, 74–92 (2013).
Article CAS PubMed PubMed Central Google Scholar
Knodel, M. M. et al. Synaptic bouton properties are tuned to best fit the prevailing firing pattern. Front. Comput. Neurosci. 8, 101 (2014).
Article PubMed PubMed Central Google Scholar
Rothman, J. S., Kocsis, L., Herzog, E., Nusser, Z. & Silver, R. A. Physical determinants of vesicle mobility and supply at a central synapse. Elife 5, 15133 (2016).
Article Google Scholar
Dittrich, M. et al. An excess-calcium-binding-site model predicts neurotransmitter release at the neuromuscular junction. Biophys. J. 104, 2751–2763 (2013).
Article CAS PubMed PubMed Central Google Scholar
Haspinger, D. C., Klinge, S. & Holzapfel, G. A. Numerical analysis of the impact of cytoskeletal actin filament density alterations onto the diffusive vesicle-mediated cell transport. PLoS Comput. Biol. 17, 1008784 (2021).
Article Google Scholar
Kim, T. & Tanaka-Yamamoto, K. Postsynaptic stability and variability described by a stochastic model of endosomal trafficking. Front. Cell. Neurosci. 13, 72 (2019).
Article CAS PubMed PubMed Central Google Scholar
Becherer, U. & Rettig, J. Vesicle pools, docking, priming, and release. Cell. Tissue Res. 326, 393–407 (2006).
Article PubMed Google Scholar
Rizo, J., Sari, L., Qi, Y., Im, W. & Lin, M. M. All-atom molecular dynamics simulations of synaptotagmin-snare-complexin complexes bridging a vesicle and a flat lipid bilayer. eLife 11, 76356 (2022).
Article Google Scholar
Marrink, S. J. et al. Computational modeling of realistic cell membranes. Chem. Rev. 119, 6184–6226 (2019).
Article CAS PubMed PubMed Central Google Scholar
Sharma, S. & Lindau, M. Molecular mechanism of fusion pore formation driven by the neuronal snare complex. Proc. Natl Acad. Sci. 115, 12751–12756 (2018).
Article CAS PubMed PubMed Central Google Scholar
Soykan, T. et al. Synaptic vesicle endocytosis occurs on multiple timescales and is mediated by formin-dependent actin assembly. Neuron 93, 854–866 (2017).
Article CAS PubMed Google Scholar
Südhof, T. C. The synaptic vesicle cycle: a cascade of protein–protein interactions. Nature 375, 645–653 (1995).
Article PubMed Google Scholar
Gillespie, D. T. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81, 2340–2361 (1977).
Article CAS Google Scholar
Hepburn, I., Chen, W., Wils, S. & De Schutter, E. STEPS: efficient simulation of stochastic reaction–diffusion models in realistic morphologies. BMC Syst. Biol. 6, 36 (2012).
Article PubMed PubMed Central Google Scholar
Hepburn, I., Cannon, R. & De Schutter, E. Efficient calculation of the quasi-static electrical potential on a tetrahedral mesh and its implementation in steps. Front. Comput. Neurosci. 7, 129 (2013).
PubMed PubMed Central Google Scholar
Chen, W. & De Schutter, E. Parallel steps: large scale stochastic spatial reaction-diffusion simulation with high performance computers. Front. Neuroinformatics 11, 13 (2017).
Article Google Scholar
Hepburn, I., Chen, W. & De Schutter, E. Accurate reaction-diffusion operator splitting on tetrahedral meshes for parallel stochastic molecular simulations. J. Chem. Phys. 145, 054118 (2016).
Article CAS PubMed Google Scholar
Isaacson, S. A. The reaction-diffusion master equation as an asymptotic approximation of diffusion to a small target. SIAM J. Appl. Math. 70, 77–111 (2009).
Article Google Scholar
Antunes, G. & De Schutter, E. A stochastic signaling network mediates the probabilistic induction of cerebellar long-term depression. J. Neurosci. 32, 9288–9300 (2012).
Article CAS PubMed PubMed Central Google Scholar
Hepburn, I. et al. A model of induction of cerebellar long-term depression including rkip inactivation of raf and mek. Front. Mol. Neurosci. 10, 19 (2017).
Article PubMed PubMed Central Google Scholar
Zamora Chimal, C. G. & De Schutter, E. Ca2+ requirements for long-term depression are frequency sensitive in purkinje cells. Front. Mol. Neurosci. 11, 438 (2018).
Article CAS PubMed PubMed Central Google Scholar
Matolcsi, M. & Giordano, N. A novel explanation for observed camkii dynamics in dendritic spines with added egta or bapta. Biophys. J. 108, 975–985 (2015).
Article CAS PubMed PubMed Central Google Scholar
Anwar, H., Hepburn, I., Nedelescu, H., Chen, W. & De Schutter, E. Stochastic calcium mechanisms cause dendritic calcium spike variability. J. Neurosci. 33, 15848–15867 (2013).
Article CAS PubMed PubMed Central Google Scholar
Anwar, H. et al. Dendritic diameters affect the spatial variability of intracellular calcium dynamics in computer models. Front. Cell. Neurosci. 8, 168 (2014).
Article PubMed PubMed Central Google Scholar
Mohapatra, N. et al. Spines slow down dendritic chloride diffusion and affect short-term ionic plasticity of gabaergic inhibition. Sci. Rep. 6, 1–14 (2016).
Article Google Scholar
Schelker, M. et al. Viral rna degradation and diffusion act as a bottleneck for the influenza a virus infection efficiency. PLoS Comput. Biol. 12, 1005075 (2016).
Article Google Scholar
Denizot, A., Arizono, M., Nägerl, U. V., Soula, H. & Berry, H. Simulation of calcium signaling in fine astrocytic processes: Effect of spatial properties on spontaneous activity. PLoS Comput. Biol. 15, 1006795 (2019).
Article Google Scholar
Stillman, N. R. et al. Evolutionary computational platform for the automatic discovery of nanocarriers for cancer treatment. NPJ Comput. Mater. 7, 1–12 (2021).
Article Google Scholar
Chen, W., Hepburn, I., Martyushev, A., De Schutter, E. Modeling neurons in 3d at the nanoscale. In Computational Modelling of the Brain: Modelling Approaches to Cells, Circuits and Networks, 3–24 (Springer, Cham, 2022).
Gallimore, A.R., Hepburn, I., Rizzoli, S., Schutter, E.D. Dynamic regulation of vesicle pools in a detailed spatial model of the complete synaptic vesicle cycle. bioRxiv https://arxiv.org/abs/https://www.biorxiv.org/content/early/2023/08/06/2023.08.03.551909.full.pdf (2023).
Park, C., Jung, S. & Park, H. Single vesicle tracking for studying synaptic vesicle dynamics in small central synapses. Curr. Opin. Neurobiol. 76, 102596 (2022).
Article CAS PubMed Google Scholar
Milovanovic, D. & De Camilli, P. Synaptic vesicle clusters at synapses: a distinct liquid phase? Neuron 93, 995–1002 (2017).
Article CAS PubMed PubMed Central Google Scholar
Rizzoli, S. O. & Betz, W. J. Synaptic vesicle pools. Nat. Rev. Neurosci. 6, 57–69 (2005).
Article CAS PubMed Google Scholar
Wilhelm, B. G. et al. Composition of isolated synaptic boutons reveals the amounts of vesicle trafficking proteins. Science 344, 1023–1028 (2014).
Article CAS PubMed Google Scholar
Hirokawa, N. & Takemura, R. Molecular motors and mechanisms of directional transport in neurons. Nat. Rev. Neurosci. 6, 201–214 (2005).
Article CAS PubMed Google Scholar
Chenouard, N., Xuan, F. & Tsien, R. W. Synaptic vesicle traffic is supported by transient actin filaments and regulated by pka and no. Nat. Commun. 11, 1–16 (2020).
Article Google Scholar
Wang, Z. et al. Myosin vb mobilizes recycling endosomes and ampa receptors for postsynaptic plasticity. Cell 135, 535–548 (2008).
Article CAS PubMed PubMed Central Google Scholar
Vale, R. D. The molecular motor toolbox for intracellular transport. Cell 112, 467–480 (2003).
Article CAS PubMed Google Scholar
Kodera, N., Yamamoto, D., Ishikawa, R. & Ando, T. Video imaging of walking myosin v by high-speed atomic force microscopy. Nature 468, 72–76 (2010).
Article CAS PubMed Google Scholar
Rief, M. et al. Myosin-v stepping kinetics: a molecular model for processivity. Proc. Natl Acad. Sci. 97, 9482–9486 (2000).
Article CAS PubMed PubMed Central Google Scholar
Mori, T., Vale, R. D. & Tomishige, M. How kinesin waits between steps. Nature 450, 750–754 (2007).
Article CAS PubMed Google Scholar
Stepp, W.L. & Ökten, Z. Resolving kinesin stepping: one head at a time. Life Sci. Alliance 2, 201900456 (2019).
Schnitzer, M. J. & Block, S. M. Kinesin hydrolyses one atp per 8-nm step. Nature 388, 386–390 (1997).
Article CAS PubMed Google Scholar
Veigel, C., Wang, F., Bartoo, M. L., Sellers, J. R. & Molloy, J. E. The gated gait of the processive molecular motor, myosin v. Nat. Cell. Biol. 4, 59–65 (2002).
Article CAS PubMed Google Scholar
Hilfiker, S. et al. Synapsins as regulators of neurotransmitter release. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 354, 269–279 (1999).
Article CAS Google Scholar
Rizzoli, S. O. Synaptic vesicle recycling: steps and principles. EMBO J. 33, 788–822 (2014).
Article CAS PubMed PubMed Central Google Scholar
Ghosh, A., Samuel, J. & Sinha, S. A “gaussian" for diffusion on the sphere. EPL Europhys. Lett. 98, 30003 (2012).
Article Google Scholar
Südhof, T. C. Neurotransmitter release: the last millisecond in the life of a synaptic vesicle. Neuron 80, 675–690 (2013).
Article PubMed Google Scholar
Fesce, R., Grohovaz, F., Valtorta, F. & Meldolesi, J. Neurotransmitter release: fusion or ‘kiss-and-run’? Trends Cell. Biol. 4, 1–4 (1994).
Article CAS PubMed Google Scholar
Harata, N. C., Aravanis, A. M. & Tsien, R. W. Kiss-and-run and full-collapse fusion as modes of exo-endocytosis in neurosecretion. J. Neurochem. 97, 1546–1570 (2006).
Article CAS PubMed Google Scholar
Alabi, A. A. & Tsien, R. W. Perspectives on kiss-and-run: Role in exocytosis, endocytosis, and neurotransmission. Ann. Rev. Physiol. 75, 393–422 (2013).
Article CAS Google Scholar
Ren, L. et al. The evidence for open and closed exocytosis as the primary release mechanism. Q. Rev. Biophys. 49, 12 (2016).
Article Google Scholar
Li, X., Majdi, S., Dunevall, J., Fathali, H. & Ewing, A. G. Quantitative measurement of transmitters in individual vesicles in the cytoplasm of single cells with nanotip electrodes. Angew. Chem. Int. Ed. 54, 11978–11982 (2015).
Article CAS Google Scholar
Bonifacino, J. S. & Glick, B. S. The mechanisms of vesicle budding and fusion. Cell 116, 153–166 (2004).
Article CAS PubMed Google Scholar
Kaksonen, M. & Roux, A. Mechanisms of clathrin-mediated endocytosis. Nat. Rev. Mol. Cell. Biol. 19, 313–326 (2018).
Article CAS PubMed Google Scholar
Südhof, T. C. The synaptic vesicle cycle. Ann. Rev. Neurosci. 27, 509–547 (2004).
Article PubMed Google Scholar
Rettig, J. & Neher, E. Emerging roles of presynaptic proteins in ca++-triggered exocytosis. Science 298, 781–785 (2002).
Article CAS PubMed Google Scholar
Denker, A. & Rizzoli, S. Synaptic vesicle pools: An update. Front. Synaptic Neurosci. 2, 135 (2010).
PubMed PubMed Central Google Scholar
Kamin, D. et al. High- and low-mobility stages in the synaptic vesicle cycle. Biophys. J. 99, 675–684 (2010).
Article CAS PubMed PubMed Central Google Scholar
Hirokawa, N., Sobue, K., Kanda, K., Harada, A. & Yorifuji, H. The cytoskeletal architecture of the presynaptic terminal and molecular structure of synapsin 1. J. Cell. Biol. 108, 111–126 (1989).
Article CAS PubMed Google Scholar
Lingwood, D. & Simons, K. Lipid rafts as a membrane-organizing principle. Science 327, 46–50 (2010).
Article CAS PubMed Google Scholar
Pike, L. J. Rafts defined: a report on the keystone symposium on lipid rafts and cell function. J. Lipid Res. 47, 1597–1598 (2006).
Article CAS PubMed Google Scholar
Hernández-Adame, P. L. et al. Determination of the size of lipid rafts studied through single-molecule fret simulations. Biophys. J. 120, 2287–2295 (2021).
Article PubMed PubMed Central Google Scholar
Simons, K. & Toomre, D. Lipid rafts and signal transduction. Nat. Rev. Mol. Cell. Biol. 1, 31–39 (2000).
Article CAS PubMed Google Scholar
Edwards, R. H. Neurotransmitter release: variations on a theme. Curr. Biol. 8, 883–885 (1998).
Article Google Scholar
Strobl, S., Formella, A. & Pöschel, T. Exact calculation of the overlap volume of spheres and mesh elements. J. Comput. Phys. 311, 158–172 (2016).
Article Google Scholar
Hart, P. E., Nilsson, N. J. & Raphael, B. A formal basis for the heuristic determination of minimum cost paths. IEEE Trans. Syst. Sci. Cybern. 4, 100–107 (1968).
Article Google Scholar
Isaacson, S. A. & Peskin, C. S. Incorporating diffusion in complex geometries into stochastic chemical kinetics simulations. SIAM J. Sci. Comput. 28, 47–74 (2006).
Article Google Scholar
Ahrens, J., Geveci, B., Law, C. ParaView: An End-User Tool for Large Data Visualization. In Visualization Handbook, 717–731. Elsevier Inc. https://www.sciencedirect.com/book/9780123875822/visualization-handbook (2005).
Hepburn, I., Lallouette, J. Vesicle and reaction-diffusion hybrid modeling with STEPS [Data set]. Zenodo, https://doi.org/10.5281/zenodo.10901994 (2024).
Hepburn, I., Lallouette, J., Chen, W. Cns-oist/steps: Steps 5.0.1. Zenodo, https://doi.org/10.5281/zenodo.10910455 (2024).

Download references

Acknowledgements

Research reported in this publication was supported by the Okinawa Institute of Science and Technology Graduate University (OIST), and by Kakenhi Grant-in-Aid for Scientific Research Project Number 20K06593. We thank Tristan Carel of the Blue Brain Project, École Polytechnique Fédérale de Lausanne for his contributions to the vesicle code.

Author information

Authors and Affiliations

Computational Neuroscience Unit, Okinawa Institute of Science and Technology, 1919-1 Tancha, Onna-son, Okinawa, Japan
Iain Hepburn, Jules Lallouette, Weiliang Chen, Andrew R. Gallimore, Sarah Y. Nagasawa-Soeda & Erik De Schutter

Authors

Iain Hepburn
View author publications
You can also search for this author in PubMed Google Scholar
Jules Lallouette
View author publications
You can also search for this author in PubMed Google Scholar
Weiliang Chen
View author publications
You can also search for this author in PubMed Google Scholar
Andrew R. Gallimore
View author publications
You can also search for this author in PubMed Google Scholar
Sarah Y. Nagasawa-Soeda
View author publications
You can also search for this author in PubMed Google Scholar
Erik De Schutter
View author publications
You can also search for this author in PubMed Google Scholar

Contributions

E.D.S. conceptualized and led the study. I.H. led software development. J.L. and W.C. contributed to software development and I.H., J.L. and W.C. did pre-release testing and debugging. I.H. and J.L. wrote the Python model scripts, ran simulations and collected data. J.L. wrote the Blender visualization module. A.R.G. and S.Y.N.S. designed features of the software and tested the implementations. I.H., J.L., W.C. and E.D.S. wrote the manuscript. All authors read and approved the final version.

Corresponding author

Correspondence to Erik De Schutter.

Ethics declarations

Competing interests

The authors declare no competing interests.

Peer review

Peer review information

Communications Biology thanks Xianchan Li, Olivia Eriksson, Andrei Kramer and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Johannes Stortz. A peer review file is available.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Peer review file

Supplementary Information

Description of additional supplementary files

Supplementary Movie 1

Supplementary Movie 2

Supplementary Movie 3

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Cite this article

Hepburn, I., Lallouette, J., Chen, W. et al. Vesicle and reaction-diffusion hybrid modeling with STEPS. Commun Biol 7, 573 (2024). https://doi.org/10.1038/s42003-024-06276-5

Download citation

Received: 11 September 2023
Accepted: 01 May 2024
Published: 15 May 2024
DOI: https://doi.org/10.1038/s42003-024-06276-5