Reconstruction of endosomal organization and function by a combination of ODE and agent-based modeling strategies

Mayorga, Luis S.; Cebrian, Ignacio; Verma, Meghna; Hoops, Stefan; Bassaganya-Riera, Josep

doi:10.1186/s13062-018-0227-4

Research
Open access
Published: 23 November 2018

Reconstruction of endosomal organization and function by a combination of ODE and agent-based modeling strategies

Luis S. Mayorga ORCID: orcid.org/0000-0002-5995-0671¹,
Ignacio Cebrian¹,
Meghna Verma^2,3,
Stefan Hoops^2,4 &
…
Josep Bassaganya-Riera²

Biology Direct volume 13, Article number: 25 (2018) Cite this article

3904 Accesses
6 Citations
1 Altmetric
Metrics details

Abstract

Background

Reproducing cell processes using an in silico system is an essential tool for understanding the underlying mechanisms and emergent properties of this extraordinary complex biological machine. However, computational models are seldom applied in the field of intracellular trafficking. In a cell, numerous molecular interactions occur on the surface or in the interior of membrane-bound compartments that continually change position and undergo dynamic processes of fusion and fission. At present, the available simulation tools are not suitable to develop models that incorporate the dynamic evolution of the cell organelles.

Results

We developed a modeling platform combining Repast (Agent-Based Modeling, ABM) and COPASI (Differential Equations, ODE) that can be used to reproduce complex networks of molecular interactions. These interactions occur in dynamic cell organelles that change position and composition over the course of time. These two modeling strategies are fundamentally different and comprise of complementary capabilities. The ODEs can easily model the networks of molecular interactions, signaling cascades, and complex metabolic reactions. On the other hand, ABM software is especially suited to simulate the movement, interaction, fusion, and fission of dynamic organelles. We used the combined ABM-ODE platform to simulate the transport of soluble and membrane-associated cargoes that move along an endocytic route composed of early, sorting, recycling and late endosomes. We showed that complex processes that strongly depend on transport can be modeled. As an example, the hydrolysis of a GM2-like glycolipid was programmed by adding a trans-Golgi network compartment, lysosomal enzyme trafficking, endosomal acidification, and cholesterol processing to the simulation model.

Conclusions

The model captures the highly dynamic nature of cell compartments that fuse and divide, creating different conditions for each organelle. We expect that this modeling strategy will be useful to understand the logic underlying the organization and function of the endomembrane system.

Reviewers

This article was reviewed by Drs. Rafael Fernández-Chacón, James Faeder, and Thomas Simmen.

Background

Depict the fact that we have a detailed catalog of factors and processes necessary for intracellular transport, our understanding of the underlying molecular mechanism is still fragmentary and qualitative. Generally, models (following Gunawardena’s definition for models: “some form of symbolic representation of our assumptions about reality” [1]) related to the movement of compounds along the endocytic pathway are informal (“one in which the symbols are mental, verbal, or pictorial, perhaps a scrawl of blobs and arrows on the whiteboard”; see Fig. 1a for an example). These informal models are powerful tools to envision potential molecular mechanisms, general cell strategies, and qualitative behaviors. However, they are difficult to validate (or falsify) because they are so vague that several alternative interpretations exists. For example, the brown arrows in Fig. 1a symbolizes the recycling and degradative pathways of an internalized molecule; however, they can also be interpreted as time and/or space fluxes. Moreover, it is not clear whether the depicted compartments remain stable or evolve as the cargoes are traveling along these two pathways. In brief, it is challenging to propose or construct models that describe the flux of cargoes through dynamic structures that continuously change shape, position, and composition.

Mathematical models have proven to be useful tools for understanding the underlying mechanisms and emergent properties in cell biology. Excellent tools are available to represent and model complex biological pathways [2, 3]. However, in the cell, many molecular interactions occur on the surface or the interior of membrane-bound compartments that continuously change position and undergo the dynamic processes of fusion and fission. At present, the available simulation tools fail to incorporate the dynamic evolution of organelles. In particular, it is troublesome to simulate the trafficking of macromolecules along the endocytic and secretory pathways [4]. In this scenario, our aim is to develop a platform for intracellular transport modeling that combines Agent-Based Modeling (ABM) and Ordinary Differential Equations (ODE) that are able to handle the dynamic nature of organelles along with the complex network of molecular interactions that occur in individual structures. In particular, we focused on the endocytic route.

Endocytosis has a central role in the physiology of the cell. It is required for the uptake of macromolecules to be processed in lysosomes, the internalization and recycling of receptors involved in several cell signaling cascades, and a myriad of diverse functions that occur in these membrane-bound organelles. By functional and morphological criteria, the endocytic route is divided into several compartments. Early endosomes receive the material endocytosed by different mechanisms (clathrin-dependent and independent endocytosis, macropinocytosis, phagocytosis, to name a few). From this early compartment, the compounds can directly recycle back to the plasma membrane (early recycling) or can be transported to the sorting endosomes where the molecules are directed to their specific intracellular destinations. Several compounds are directed to the endosomal recycling center, located near the nucleus, and responsible for the slow recycling to the cell surface. If the endocytosed materials are not sorted out from early endosomes, their final destination is the late endosomal/lysosomal compartment triggered by the maturation of early to late endosomes [5]. The presence of hydrolytic enzymes and acidic environment in the lysosomes promote the digestion of most of the biological macromolecules.

To develop a formal model of the endocytic pathway, we utilized the concept that endosomes are characterized by their limiting membranes. Membrane domains are relatively stable structures that do not necessarily mix when residing side by side in the same organelle [6, 7]. Membrane domains can thus be considered as the key building blocks of the endomembrane system of the cell. Transport can be envisioned as the change in the limiting membrane of organelles carrying a cargo by means of mixing (fusion of organelles) and sorting (budding and fission of vesicles and tubules) processes. In addition, complex networks of molecular interactions change the protein and lipid composition of the limiting membrane of the organelles [8, 9]. Soluble cargoes that do not interact with the membranes are transported in the lumen of the organelles. During the budding of vesicles and tubules, these cargoes are partitioned in volumes proportional to those of the splitting structures; hence, very little is incorporated in the newly-formed membrane-rich structures whereas a large proportion is retained in the remaining round organelle that carries most of the total volume. Since tubules are key components of the recycling pathways, most soluble cargoes are delivered to the lysosomes [10]. On the contrary, membrane-bound factors travel through the endomembrane system according to their affinity for membrane domains that are sorted out during the budding of vesicles or tubules. Membrane association depends on specific interactions of tags coded in the cargo molecule with one or more of the large set of adaptor proteins present in the cell [11, 12]. In addition, lipids on the membranes undergo lateral phase separations, and membrane-anchored factors can also be recruited to specific lipids domains [13]. Depending on their protein and lipid composition, different membrane domains bind membrane-deforming protein complexes, such as COPs, clathrin, sorting nexins, ESCORT, etc. [14,15,16,17]. The deformations lead to the sorting of membrane domains by budding of tubules or vesicles, or the isolation in intraluminal vesicles that are now separated from the membrane of the original organelle, carrying a subset of cargo molecules. Membrane domains depend on the coordinated activity of multiple factors, including several members of the Rab family [7, 18]. Perturbation of Rab functions, dramatically alters the intracellular transport of most of the cargoes [19, 20]. Rabs are small GTPases that act as molecular “on/off” switches, cycling between inactive (GDP-bound) and active (GTP-bound) states. They can be organized in cascades adding flexibility and regulation to the identity of membrane domains [5, 6, 21]. Furthermore, Rab proteins localize to specific membrane-bound compartments; hence, they can be used to identify specific membrane domains.

The simulation implemented using the modeling platform reproduced many aspects of the endocytic transport. Enzymes, molecular pumps, structural proteins, and lipids can be included in the simulation to study complex signaling and metabolic networks. The model captures the highly dynamic nature of cellular structures that fuse and divide, creating particular conditions for molecular interactions that are different for each organelle. We expect that this combination of modeling strategies will be useful for cell biologists to generate formal models that can be used to unveil the essential mechanisms underlying the highly efficient transport processes observed in living cells.

Results

Brief model description

In this section, we provide a simplified and intuitive description of the model. More details can be found in Methods and in the notes included in the Repast and COPASI code that are freely accessible at the Git repository https://github.com/ihem-institute/immunity/tree/LipidMetabolism.

Organelles carrying four different membrane domains were modeled in Repast (ABM platform). For simplicity, we can think of organelles limited by RabA as early endosomes, by RabB as sorting endosomes, by RabC as recycling endosomes, and by RabD as late endosomes/lysosomes (Fig. 1a). The organelles, carrying membrane and soluble content were allowed to move, interact and change their shape with time. To setup the simulation, the Repast Builder generates a set of initial organelles according to the parameters and characteristics selected by the user (Fig. 1d and Additional file 1). Repast then runs a loop of actions for all the agents in the model as listed in Fig. 1d and described in Methods. The user can control the model by updating parameters or pausing the simulation. The model continuously generates records of the evolution of the system as requested by the user (example in Additional files 2 and 3).

Shape and movement

Endosomes in the model have a membrane area that embraces a luminal volume and are represented by ellipsoids (spheres for round structures and elongated ellipsoids for non-round endosomes, depending on the volume/area ratio; Figs. 1b, 2, and movie in Additional file 4). The size of the organelles is represented in a nm scale. The endosomes are allowed to move randomly except when near a microtubule (MT, light blue straight lines in the model, Figs. 1a, 2, and movie in Additional file 4). Round and tubular organelles on MT, move to the surface or to the center of the cell according to their membrane domains at a speed of 1 μm/sec.

Fusion, fission, and luminal vesicles

Two closely spaced endosomes with compatible membrane domains are allowed to fuse and form a single organelle. This organelle carries all the membrane and soluble components of the original endosomes. The endosomal compatibilities are listed in Table 1. Endosomes with enough membrane can also split, and form tubules (40 nm diameter) carrying a single membrane domain. The soluble content is distributed proportionally to the volume and the membrane components proportionally to the area of the newly formed structures. The membrane-associated cargoes may have affinities for some specific membrane domains, and will partition during fission according to these affinities. Endosomes can also form internal vesicles, thereby decreasing their limiting area and increasing their volume (the number of internal vesicles is shown as a number inside each endosome, Figs. 1b and 2). Note that the area, volume, Rab, membrane, and soluble contents were preserved during fusion and fission events.

Table 1 Parameters used to characterize the organelles simulated in the model

Full size table

New endosomes

When new RabA organelles were formed, they were modeled as an endocytic event, budding from the cell surface and carrying soluble and membrane-bound components.

Recycling

When endosomes carrying RabA or RabC membrane domains were proximal to the plasma membrane, they could release their membrane and soluble content (recycled material), that was incorporated into the plasma membrane or released to the medium.

Lysosomal digestion

Endosomes carrying RabD domains digest part of their content and membrane.

Rab dynamics

We generated a set of ODE in COPASI to model the self-preserving behavior of membrane domains and a Rab switch homologous to the Rab5/Rab7 conversion to regulate the transformation of early (RabA) to late (RabD) endosomes (Fig. 1c, Table 3, and Additional file 5) [22]. COPASI and Repast communication was achieved using a previously published method [23]. Hence, in the simulations, Rab dynamics was controlled by differential equations (COPASI) and the rest of the model by ABM (Repast).

A global analysis of the endocytic route proposed shows four stable compartments that communicate with each other by fusion and fission events (according to membrane domain compatibilities shown in Table 1. Notice that homotypic fusion has the highest probability for all compartments. The exchange of material during fusion/fission processed is in principle bidirectional. However, the flux can be orientated according to how the molecules distribute during fission events. On top of this organization, two irreversible steps strongly influence transport: the maturation of early to late endosomes and the recycling to the plasma membrane. In the present model, only the content is recycled; hence recycling endosome domains are preserved. In contrast, maturation caused the consumption of early endosomes and the building up of late endosomes. To compensate for this, new early endosomes are constantly formed and lysosomal digestion consumes RabD domains.

Model characterization

In ABM models, during a tick interval, all agents are interrogated about performing or not a set of actions (Fig. 1d and Methods). The tick duration was calibrated with the fastest process in the model (movement of organelles on MT; 1 tick = 0.06 s). The model was stable for more than 180.000 ABM ticks corresponding to about 3 h in “cellular” time. The characteristics of the resulting organelles correspond to a rather standard endocytic pathway with organelles moving in and out the “cell” surface by random movements or along microtubules (Fig. 2 and movie in Additional file 4). Large round-shaped structures accumulated opposite to the plasma membrane (in the nucleus side). Fusion and fission of the organelles triggered the mixing and segregation of membrane domains, while the RabA/RabD switch guaranteed the eventual maturation of early to late structures. Events of fusion, fission, maturation and cargo transfer are shown in the inset of Fig. 2a-c and explained in the figure’s legend.

If a single endosome was generated, a RabA organelle matured to RabD according to the molecular switch modeled by differential equations in COPASI (Fig. 1c and Table 3). With time, RabD was incorporated into the organelle and finally the switch occurred (Fig. 3a). When other Rabs were also present, the process was more complex and the prevailing Rab tended to displace less abundant Rabs. Fig. 3b shows what happened with an organelle that initially had 45% of its surface covered with RabB and 55% with RabA. In this condition, RabB disappeared and RabA covered the endosome to, later on, switched to RabD. However, notice that the switch kinetics was altered by the presence of RabB. On the contrary, if 55% of the organelle was initially covered with RabB, the endosome rapidly evolved to an all-RabB structure (Fig. 3c).

The scenario was far more complex when other endosomes were included in the model. The molecular switch was profoundly affected when a new membrane domain was incorporated by fusion with other endosomes or when fission segregated a complete membrane domain from the organelle. Examples of RabA/RabD switches are shown in Fig. 3 d to g. In Fig. 3d, a fusion occurred just after the switch, causing a large increase in the endosome size (size are shown as relative to the initial endosome radius underneath the corresponding switch process). In Fig. 3e, there was a fusion (size increase) and fission (size decreased) after which, the endosome matured to a RabD structure. In Fig. 3f, there were several fusion and fission events that profoundly affected the switch kinetics. Finally, in Fig. 3g, the endosome was first converted to a RabB structure and then to a RabD organelle by fusion and fission processes.

The behavior of this complex scenario can be analyzed in different ways. For example, although in Fig. 3a to c it is clear that a single endosome was analyzed, it is not obvious which endosome is followed in Fig. 3d to g. This is due to the occurrence of fission events that make it unclear which structure is described. To avoid this ambiguity, markers, modeled as soluble or membrane-associated cargoes, were loaded in the endosomes. The analysis of the organelles that contain a specific marker during the course of the simulation provided meaningful information about the transport events for a specific cargo. In Fig. 3d to g, a membrane marker was included and the Rab composition of the endosome carrying the marker was followed along the time. Fig. 3h and i, show the characteristics of organelles carrying a membrane and a soluble cargo along the course of the 30 min simulation. Both markers were initially loaded in the same RabA structure. Notice that both sensed the same initial steps of a RabA/RabD switch, but before maturation was completed, fusion and fission events transformed the organelle in a RabB structure, from where the fate of the two markers diverged. The membrane marker was segregated in a RabC organelle and for a while kept moving between endosomes carrying RabB and RabC domains, until it was recycled to the plasma membrane from a RabC endosome. Notice the changes in the organelle size that oscillated between 40 and 100 nm radius (Fig. 3j). The soluble marker instead, was transported from a RabB endosome to a large RabD organelle that kept growing by fusion with other RabD structures (Fig. 3j).

A different way to analyze the simulation is to take all the organelles at a given time and classify the structures according to the prevailing Rab domain as early, sorting, recycling or late endosomes and to analyze the characteristics of these structures. Histograms in Fig. 4a and b show the size and internal vesicle content distribution of all organelles classified according to their prevailing Rab domain after 30 min of simulation. Notice that all the organelles have some small endosomes (radius in the 20–40 nm range); however, large structures (radius > 300 nm) were more frequently observed in late endosomes (RabD) and were absent in the early endosomes (RabA). The size of the largest endosomes carrying RabB or RabC domains had an intermedia radius in the range of 100–200 nm (Fig. 4a). Internal vesicles were observed in all endosome types, but they were more frequently seen in RabB and RabD structures (Fig. 4b).

Another analysis that we performed includes the calculation of the amount of a cargo associated with the different compartments at each time point of the simulation (see Methods). For example, the total number of internal vesicles in the simulation can be classified as being associated to either RabA, RabB, RabC or RabD structures (Fig. 4c). Notice that most of the vesicles accumulated in the RabD endosomes. About 500 internal vesicles were observed in the RabD structures at time 30 min, whereas, RabB and RabC had about 100 vesicles (Fig. 4c, left panel) and RabA had only 1 vesicle (Fig. 4c, right panel). In the expanded scale shown at the right, it can be seen that some internal vesicles recycled to the plasma membrane, a process that mimics the release of exosomes from multivesicular bodies (recycled, Fig. 4c, right panel).

Transport of soluble and membrane-bound cargoes

In the complex scenario simulated, transport of a cargo molecule is profoundly affected by its behavior during fission events. Like proteins carrying a specific tag in their sequence, cargo moved in this artificial endocytic route depending on their affinity for membrane domains that conditioned their segregation in one or other structure during the budding and fission of tubules. On top of these processes, the maturation of early endosomes will transport to late endosomes all cargoes that are not sorted out to different organelles. Soluble cargoes are not recognized by membranes and accumulated in the vesicular (large volume) structure that remains after the formation of tubules. Moreover, large cargoes (larger than 40 nm diameter, like intraluminal vesicles) are excluded from tubules. These cargoes are prompted to be finally delivered to large late endosomes (Fig. 4c).

To analyze the behavior of different cargoes, a soluble (dextran-like) and two membrane-bound cargoes were loaded into early endosomes (RabA) and the model was allowed to evolve for 30 min. The two membrane proteins simulated included: i) a transferrin-like cargo (Tf) that could recycle from early endosomes and recycling endosomes, and ii) a class I major histocompatibility complex -like cargo (MHC-I) that needed to be sorted to RabB structures before being recycled from RabC endosomes. Both proteins were provided with membrane affinities for RabB and RabC domains (Table 2). This means that they were sorted during fission in structures rich in these Rab domains. To analyze the movement of the cargoes, the percentage of each molecule in organelles with different membrane domains was calculated and plotted in Fig. 5. The data showed the average of 8–10 simulations and reproduced the well-characterized kinetics of pulse-chase experiments for these ligands. A clear difference was observed between the two membrane proteins Tf and MHC-I, which recycled to the plasma membrane, and dextran that accumulated in late endosomes (Fig. 5a to c). The possibility of recycling from early endosomes made the kinetics of Tf recycling more efficient than that of MHC-I. Experimental data have shown that perturbation of sorting endosomes (for example, knocking down Rab22a) profoundly affect the recycling of MHC-I molecules [24, 25] and also, although to a lesser extent, that of Tf [26]. In contrast, transport of soluble markers was only marginally altered [27]. In the simulation, when we knocked out RabB structures, these observations were reproduced fairly well. Tf recycling was strongly diminished, but not abolished whereas the recycling of MHC-I was completely blocked. Transport of dextran to late endosomes was only slightly delayed (Fig. 5d to f).

Table 2 List of molecules included in the simulation and their affinities for Rab domains

Full size table

In a different set of simulations, cytosolic RabD was depleted, preventing the maturation of RabA structures. Since RabD organelles were still present in the model, we homologate this condition to a RabD knockdown. The decreased cytosolic RabD inhibited the accumulation of dextran in late compartments and caused a dramatic increase of the fluid phase marker in sorting and recycling endosomes (Fig. 5h). In contrast, Tf recycling was not impaired. The effects observed closely reproduce what is experimentally observed when Rab7 function is perturbed [28]. The appearance of large sorting and recycling endosomes was evident (compare Fig. 5i with Fig. 2), a phenotype that has been described in cells depleted of Rab7 or expressing a dominant negative mutant of this protein [29].

In all these simulations, the endocytic markers were incorporated in RabA compartments at time zero. In a different protocol, RabA endosomes were empty and uptake was allowed in the 3–6 min time period. Under these conditions, transferrin and dextran were recycled and transported to late endosomes, respectively, with similar kinetics as observed in the time zero protocols (Additional file 6). These results indicate that transport was not strongly affected by the initial conditions of the system. In the figure shown in Additional file 6, six individual runs of the same simulation are shown. Given the stochastic nature of transport, all simulations present particular uptake events and cargo transport along time. However, in all simulations, there was a consistent recycling of transferrin and transport of dextran to late endosomes.

Exploring the power of adding molecular functionality to the simulation

The simplified endocytic route modeled here and the mechanisms for cargo transport proposed reproduced quite well the behavior of several molecules transported along this pathway. But the real strength of this simulation platform is the possibility of easily adding molecular functionality to the structures. For example, a set of differential equations was modeled in COPASI to simulate endosomal acidification. Each organelle was allowed to pump proton into the lumen by means of a membrane-associated vATPase. Leakage reactions allowing the in and out movement of protons were also included (see Methods, Table 3, and Additional file 7). The vATPase was given low affinity for RabB and RabC domains (0.0), high affinity for RabD structures (1.0) and a moderate affinity for RabA (0.7) (Table 2). With these values, the vATPase was enriched in the early/late endosome path, but some activity leaked to other compartments, a distribution that fairly mimicked what is observed in most cells [30].

Table 3 Set of reactions and kinetic functions programmed in COPASI

Full size table

After a few minutes of simulation, a pH gradient was established in the system. The analysis of this parameter in all endosomes after 30 min of simulation indicated that RabA endosomes had a pH average of 5.7 with a very broad distribution, whereas RabB and RabC were only lightly acidic with pH average of 6.3 and 6.4, respectively (Fig. 6a). On the other hand, RabD organelles had an acidic pH with an average of 5 (Fig. 6a). Following the pH sensed by membrane-associated (MHC-I-like) markers, we observed that as they moved along different organelles, on average, they were exposed to near neutral pH (6.4 ± 0.22, mean ± SD, taking into account the time remaining at each pH, N = 10, Fig. 6b). However, the probes transiently visit organelles with low pH (Fig. 6b). On the other hand, soluble markers moved rapidly to acidic organelles rich in vATPase with very few oscillations in the pH (pH ~ 5, Fig. 6b).

Endosomal pH is an important property of these organelles and plays a central role in several aspects of the endocytic route. For example, the consensus is that acid hydrolases are delivered to lysosomes by vesicular structures coming from the trans-Golgi network (TGN) that fuse with early/late endosomes [31]. The enzymes, bound to transport proteins (e.g., mannose 6-phosphate receptor, M6PR), are released into the lumen at the acidic pH of endosomes. To simulate this transport in the model, a new organelle carrying a different membrane domain was included (Fig. 1a). This organelle was given a RabE domain that could interact with RabB organelles (Table 1). A lysosomal enzyme (β-hexosaminidase A, HexA) bound to a membrane-associated transporter (M6PR) was located in the RabE structures and ODEs for membrane association and dissociation reactions were created in COPASI (Tables 2 and 3). Dissociation was promoted at pH below 6 and association was favored at higher pH. Fig 7a shows that the pH in RabE was almost neutral (6.7). Hence, the concentration of free HexA in these organelles was low (Fig. 7b). When the system was allowed to evolve for several minutes, the free enzyme accumulated almost exclusively in RabD organelles (Fig. 7b and c). The M6PR was given affinity only for RabE domains and remained during the simulation mostly in RabE structures (Fig. 7d), but in small amounts, it was found in RabB and C structures, and from there it leaked to the plasma membrane (Fig. 7e, recycled).

Another property of the endocytic route is a differential distribution of lipids. Cholesterol strongly influences the biophysical and functional properties of membranes and modulates the activity of receptors and enzymes [32]. This lipid is present in a gradient along endosomal structures. The content of cholesterol is higher in the plasma membrane and decreases from early to late endosomes. Lysosomes have a central role in the sorting of cholesterol out of the endocytic route [33, 34]. To mimic this features of the endocytic pathway, cholesterol was included as a membrane cargo in the newly formed RabA endosomes and the cholesterol of lysosomes was decreased by a reaction modeled in COPASI (Table 3). Cholesterol was given affinity for RabA (1.0) and RabB (0.5) organelles (Table 2). Under these conditions, a cholesterol gradient was spontaneously generated after a few modeling minutes (Fig. 8a).

Acid hydrolases in the lumen of late endosomes/lysosomes are responsible for the digestion of many organic molecules, including complex lipids. Defects in this process intoxicate the cell and lead to several storage diseases. Lipid metabolism is influenced by the pH, the specific set of acid hydrolases and also by the lipid composition of the membranes [35, 36]. For example, cholesterol inhibits the activity of several lipid hydrolases [36]. With the incorporation of a lysosomal enzyme (HexA), and the regulation of pH and cholesterol content, the lysosomal digestion of a complex lipid, such as the GM2 ganglioside could be simulated. With this aim, a digestion reaction, stimulated by HexA and acidic pH, and inhibited by cholesterol was programmed in COPASI (Table 3 and Additional file 7). GM2 was included in RabA endosomes as a membrane-associated cargo with affinity for RabD domains (Table 2). The continuous uptake of the glycolipid did not lead to lysosomal accumulation and a gradient was formed with recycling of undigested GM2 to the plasma membrane (Fig. 8b). Perturbation of the conditions for glycolipid digestion (alterations in pH, HexA delivery or lysosomal cholesterol content), strongly affects the efficient digestion of GM2. For example, if the lysosomal processing of cholesterol was inhibited in COPASI, RabD structures became enriched in this lipid (Fig. 8c). High levels of cholesterol in lysosomes caused a defect in GM2 metabolism and the lipid rapidly accumulated in RabD structures (Fig. 8d).

Conditions for GM2 hydrolysis were different in each endosome. The amount of lipid and hydrolytic enzyme, the pH and cholesterol content were peculiar for each organelle and changed with time. The simulation allowed to analyze by COPASI the glycolipid metabolism under the prevailing conditions of individual endosomes. Simultaneously, the Rab dynamics for the endosome can be followed along time. As examples, the lipid metabolism and Rab dynamics for four different endosomes are shown in Fig. 9. The first two endosomes acquired RabD characteristics (Fig. 9e and f), but the pH, cholesterol content and hydrolytic enzyme were not suitable for GM2 hydrolysis (Fig. 9a and b). On the contrary, the other two endosomes were rich in RabD (Fig. 9g and h) and digested the lipid but with different kinetics (Fig. 9c and d).

Discussion

We present the first formal model of a complete endocytic pathway mimicking uptake, recycling, and digestion of soluble and membrane-bound cargoes. The model can be expanded and complexity can be added to deal with other specific aspects of the pathway not included in the present version. In addition, this modeling strategy is suitable to study molecular interaction and chemical transformations occurring inside endosomes where the conditions change abruptly upon fusion and fission events. These possibilities were exemplified by adding an extra TGN-like domain in the ABM model and ODE to deal with endosomal acidification, lysosomal enzyme delivery, cholesterol processing, and glycolipid digestion. Notice that the model can be expanded to deal with additional structural and biochemical complexity. It is known that intraluminal vesicles and acidic lipids are important for GM2 digestion [35, 36]. Internal vesicles can be promoted to agents carrying specific lipids and proteins, and the network of acid hydrolases and co-activators necessary for glycolipid digestion can be included in the simulation. Also, the Golgi-endosome trafficking can be completed by adding more transport molecules for the enzymes and new membrane domains to capture the complexity of this particular transport step [12, 31]. The platform is also suitable to model processes occurring in the cytosolic face of the membranes. An elongated world as shown in Fig. 2f could be a good scenario to study molecular mechanisms of organelle movement along microtubules.

The model is also a proposal about the general principles governing endosomal recycling and lysosomal targeting. The pathway modeled requires a fine tuning between maturation, which transports in bulk soluble and membrane components to late endosomes, and fusion-fission events, able to sort out membrane factors from the maturing process and direct them to other cellular destinations, including the recycle to the cell surface. From this description, it is clear that transport requires “maturation” and “vesicular transport” steps, suggesting that both mechanisms are necessary for efficient intracellular trafficking.

In the context of fusion/fission-mediated transport, movement of membrane cargoes between endosomes having different Rab domains depends on the fusion probability between the domains and partitioning preference of the cargo when the budding of vesicles or tubules segregates the membrane domains in separate structures. The idea that endosomal sorting to other cellular destinations is necessary to avoid transport to lysosomes has been proposed as informal models several times [16, 37]. We expect that having the possibility of modeling cargoes that are directed to specific membrane domains will provide a convenient platform to put together the different sorting mechanisms already described.

The random and iterative nature of fusion/fission events generates very heterogeneous sets of endosomes carrying the same Rab domain. These differences were not stable, and changed rapidly along time as new fusion and fissions occurred. Thus we conclude that endosomes with a very specific combination of conditions that might be necessary for some important events are likely to occur, with a given probability, but will not be stable and recognizable as an independent compartment.

The model also stresses the importance of low-frequency events. Without the low but significant possibility of RabB/RabD fusion, the late endosome compartment would be an irreversible sink, without any chance of rescuing material miss-sorted to this compartment by maturation. Also, luminal vesicles, capable of being recycled to the extracellular medium (exosome release-like process), and the presence of M6PR at the plasma membrane occurred by a combination of rear events.

Other formal models for intracellular transport have been previously published, in most cases, as proof of principle of different transport mechanisms [38,39,40]. Most of the published models use ODE. A seminal example is the paper by Heinrich and Rapoport where they showed that stable non-identical compartments can be formed by a combination of SNARE-dependent fusion and coat-dependent fission [41]. ABM has been used to model the spatial and temporal dynamics of autophagic organelles [42], and the cellular network of mitochondria [43]. We have previously applied ABM to show that soluble cargoes can be efficiently transported by endosome fusion and asymmetric fission [10]. The molecular crowding effect and the influence of the cytoskeleton on intracellular transport have been studied using a more mechanistic version of ABM [44, 45]. A combination of ABM and ODE has been used to study mitophagy [46]. However, our model is the first to simulate a complete and stable endocytic route capable of handling the uptake, recycling, and digestion of membrane-bound and soluble molecules.

The modeling platform introduced has some shortcuts, for example: i) The scenario is at present a projection in two dimensions of the real 3D space. Repast has the possibility of modeling a three dimension space, but agents and rules would have to be adapted to a more complex and computationally expensive situation. ii) The modeled space is a rectangle, which is convenient to represent a section of a cell. However, some simulations may require specific shapes that will require to generate more complex borders. iii) At present, the form of organelles depends only on their volume and area, which are not enough to represent the shape of real organelles. iv) COPASI calculations assume a homogeneous distribution of the species, including membrane-associated molecules, which may be not convenient for some situations, especially for large structures. v) Agents at a given time cannot be stored to be used for new runs. It would be desirable to let the system reach a steady state and use this set of agents to start the simulations. We expect to solve this limitation in the near future. In brief, the present version of the platform provides a convenient modeling frame for intracellular transport, and its application to concrete trafficking problems will be important to show what limitations need to be addressed in new versions.

Conclusions

Modeling chemical reactions and molecular interaction occurring in the interior or the boundary of dynamic organelles constantly change position and composition is challenging. Our simulation relies on two complementary modeling strategies: ABM and ODE. Tools for modeling using ODE are freely available and are relatively user-friendly [2, 47]. Networks of molecular interactions, or signaling cascades can be easily modeled in COPASI. On the other hand, ABM software is ideal to model movement and interactions among organelles, including fusion and fission events. However, a drawback for ABM is that, because of its great flexibility, it is too generic to be used without some basic programming skills [48]. The connectivity between the two strategies is simple although it requires to be adjusted for different applications [23]. In the near future, we will work in a more user-friendly platform easily expandable to incorporate more domains and networks of proteins and lipids to model specific aspects of the endocytic pathway.

Modeling a particular intracellular transport process will require: i) to select the set of interacting organelles to be included, and define their properties (Rab domains, fusion compatibility, maturation, tubule tendency, migration on MT, recycling probability), ii) to select the set of molecules required for the process and define their affinity for Rab domains and their network of molecular interactions, and iii) to specify a starting set of organelles. For the GM2 hydrolysis simulation, organelles, molecules, and parameters are listed in Tables 1, 2, 3, and Additional file 1. Molecular interactions were programmed in COPASI (Additional files 5 and 7). In a first approach, parameters can be selected by educated guesses taking into account what is known about the process. Experimental values for specific steps or global behavior of the system can be used to adjust the model and narrow the parameter uncertainty.

We expect that our approach will contribute to the development of quantitative and formal hypotheses that will foster a better understanding of the highly complex and dynamic nature of the intracellular transport.

Methods

Agent-based model (Repast)

The freely available modeling platform Repast [48] was used to model agents and actions in an Eclipse environment (https://repast.github.io/). The code can be accessed from the Git repository https://github.com/ihem-institute/immunity/tree/LipidMetabolism

Ordinary differential equations (COPASI)

ODEs were programmed in COPASI [2] (http://copasi.org/). All COPASI files are included in the Git repository. COPASI and Repast interaction is achieved through the RabConvertion and LipidMetabolism classes. Basically, Repast sends initial concentrations present in each endosome to COPASI that generates a time series. A matrix with time series for each metabolite is sent back to Repast.

World

The space represented is a projection in 2D of a cytosol square of 4.5 × 4.5 μm. The size of the world can be changed using a scale factor (Fig. 2d and e). The shape of the world can also be changed, adjusting Repast parameters (Fig. 2f). The upper border corresponds to the plasma membrane and the lower border to the nucleus. The right and left borders form a continuous. Hence, the world shape corresponds to the surface of a cylinder.

Time

The tick duration was calibrated with the fastest process in the model (movement of organelles on MT; 1 tick = 0.06 s). The time for all other actions was adjusted changing the probability of being performed at each tick. Hence, actions (see below) occurred with the following frequencies every 0.06 s: Move, 1; Fusion, 0.001; Tethering, 0.06; Internal Vesicle, 0.012; New Endosomes, 0.001; Recycling, 0.06; Lysosomal digestion, 0.006; Fission, 0.15.

Agents

Each “Endosome” in the model has area and volume, Rab content, membrane content, and soluble content. Contents were expressed in mM units. About 20 molecules of a soluble molecule at 1 mM concentration will be present in a 20 nm radius vesicle. This value fits well with the reported range of membrane proteins in an average synaptic vesicle of 21 nm radius (2–70 molecules, [49]). Assuming that membrane molecules are in a homogeneous mix with soluble components is one of the limitations of the model.

The shape and size of endosomes correspond to a spheroid (prolate spheroid) with the area and volume of the endosome. The spheroids were represented as ellipses in the 2D representation of the world. Endosomes can perform the following actions:

Move

The endosomes move randomly according to their angular moment calculated from their shape and volume. Basically, large organelles move and change direction more slowly than small ones. Endosomes near the plasma membrane or the nucleus move always randomly. Organelles away from the borders and near a MT (light blue straight lines in the model) may move to the plus (plasma membrane) or minus (nucleus) end of the filament. Supported by experimental evidence, the speed of movement on MT was set to 1 μm/sec for all organelles, irrespective of their size, shape or membrane domains [50]. On MT, round and tubular organelles move to the surface or to the center of the cell according to their membrane domains. Tubules move to plasma membrane with a probability proportional to the area covered by RabA, B and C. Non-tubular structures move to the cell center with a probability proportional to the area covered by RabC, D and E. In summary, tubules move to the plasma membrane if they are rich in RabA, B or C; otherwise, they move randomly. Non-tubules move to the nucleus if they are rich in RabC, D and E; otherwise, they move randomly.

Fusion

Each endosome sensed all other endosomes at a distance less than its size (the radius of a sphere with the endosome’s volume). If nearby endosomes carry a compatible membrane domain, they fuse. Endosome compatibility was calculated as previously described [10]. In essence, the probability of fusion of single domain endosomes is specified in Table 1. For endosomes carrying more than one domain, the probability was adjusted according to the proportional area occupied by each Rab. Fused endosomes formed a single organelle carrying the area and volume and all the membrane and luminal (soluble) components of the original endosomes.

Tethering

As a step prior to fusion, nearby endosomes carrying compatible domains can form a cluster with a high probability of moving together. The rules to form a cluster were the same than for fusion and depended on the distance among endosomes and Rab domain compatibility.

Fission

Elongated endosomes have enough membrane to generate tubules. Long tubules could also divide into two tubules. Fission always generate tubules with a 20 nm radius carrying a single Rab domain. The Rab domain that was incorporated in the tubule was selected at random, taking into account the tubule forming probability of each Rab domain (specified in Table 1). Notice that RabD tubules grow from endosomes that have no other domains. The size of the tubule depends on the available surface of the Rab domain selected and on the endosome surface (the remaining surface should be enough to contain the volume not incorporated in the tubule). The smallest tubule formed was 40 nm long (almost a spherical vesicle). Soluble content distributed proportional to the volume of the formed structures, except luminal vesicles that could not enter into tubules. Membrane components were directed to the tubule only if they have more affinity for the Rab domain forming the tubule than for the Rab domains remaining. Since the content of membrane-associated cargoes was in area units, the budding tubules carried, at most, the cargoes corresponding to the tubule area. If the affinity for the tubule was the same than for the rest of the organelle, the cargo distributed proportionally to the area of the two new organelles. The area, volume, Rab, membrane, and soluble contents were preserved during fusion and fission events.

Internal vesicles

Endosomes carrying RabA, B or D domains could form internal vesicles (20 nm radius sphere). The area of the organelle decreased (the surface of the sphere was subtracted) and the volume increased (the volume of the sphere was added). The soluble content was not affected, but membrane content decreased proportionally to the internalized area. The internal vesicle was considered as a soluble content. The number of internal vesicles of each endosome was displayed as an integer number at the center of the organelles (Figs. 1a and 2).