A Conceptual Mathematical Model of the Dynamic Self-Organisation of Distinct Cellular Organelles

Formation, degradation and renewal of cellular organelles is a dynamic process based on permanent budding, fusion and inter-organelle traffic of vesicles. These processes include many regulatory proteins such as SNAREs, Rabs and coats. Given this complex machinery, a controversially debated issue is the definition of a minimal set of generic mechanisms necessary to enable the self-organization of organelles differing in number, size and chemical composition. We present a conceptual mathematical model of dynamic organelle formation based on interacting vesicles which carry different types of fusogenic proteins (FP) playing the role of characteristic marker proteins. Our simulations (ODEs) show that a de novo formation of non-identical organelles, each accumulating a different type of FP, requires a certain degree of disproportionation of FPs during budding. More importantly however, the fusion kinetics must indispensably exhibit positive cooperativity among these FPs, particularly for the formation of larger organelles. We compared different types of cooperativity: sequential alignment of corresponding FPs on opposite vesicle/organelles during fusion and pre-formation of FP-aggregates (equivalent, e.g., to SNARE clusters) prior to fusion described by Hill kinetics. This showed that the average organelle size in the system is much more sensitive to the disproportionation strength of FPs during budding if the vesicular transport system gets along with a fusion mechanism based on sequential alignments of FPs. Therefore, pre-formation of FP aggregates within the membranes prior to fusion introduce robustness with respect to organelle size. Our findings provide a plausible explanation for the evolution of a relatively large number of molecules to confer specificity on the fusion machinery compared to the relatively small number involved in the budding process. Moreover, we could speculate that a specific cooperativity which may be described by Hill kinetics (aggregates or Rab/SNARE complex formation) is suitable if maturation/identity switching of organelles play a role (bistability).


Vesicle size dependence of the fusion rate constant
Due to the hydration repulsion water molecules attached to the charged membranes must be removed to overcome the energy barrier ∆E ij of two fusing vesicles i and j. For the shaded spherical caps (∆O i and ∆O j ) in Fig. 9 which cover these molecules one must chose an appropriate parameter ∆ defining the necessary interaction cross section to be considered. In the simplest case one may expect ∆E ij ∝ (∆O i + ∆O j ) = 2πr i ∆ i + 2πr j ∆ j . For the two segments of ∆, named ∆ i and ∆ j (see Fig. 9), a reasonable approximation is given by ∆ i = r j /(r i + r j ) · ∆ and ∆ j = r i /(r i + r j ) · ∆ which gives a compact form of the energy barrier entering Eq. (2): playing the role of an phenomenological proportionality factor.

Hyper-linear fusion rates
The total fusion rate f ij of two vesicles i and j is composed of individual fusion rates: where v ij (m) describes the concentration of intermediates with m FP pairs already formed (see Fig. 3A). We assume rapid equilibria between successive intermediates v ij (m) and v ij (m + 1): where the factor (n−m)(n−m) describes the combinatorial multiplicity of pairing options of one additional FP pair and (m + 1) accounts for the possible unpairing processes. It can be easily shown that (1) can then be rewritten in a compact form: We assume that a successful fusion can be accomplished primarily if at least n FP pairings could be established (overcome of energy barrier), i.e. k 0 , ...k m , ..., k n−1 ≪ k n . The accelerating effect of SNARE pairing is accounted for by the association constant K* describing the stabilisation of the fusion intermediate by SNARE pairing. If no pre-formation of aggregates is assumed K * will certainly depend linearly on the concentration of total available FPs: K * ∝ A i A j : In the case of two independent FPs A and B, one retrieves where κ ′ accounts for the vesicle size dependent association constant K (repulsive hydration force). A specific assumption of this dependency yields Eq. (5).
In the scenario of pre-forming aggregates (see Fig. 3B) a concerted transition from v ij (0) to v ij (n) is assumed to proceed by a single step. In this case K * depends on the number of aggregates instead of single FPs giving rise to a Hill type cooperativity as described in Eq. (6).

De novo generation of unit size vesicles
The enrichment of unit vesicles with either A or B type FPs results in coatA or coatB unit vesicles, respectively. For the initial de novo generation of such vesicles as well as for the budding process from organelles (see next section) we need a quantifier to describe the degree of biased loading during these two processes. We therefore introduce an enrichment factor η > 1 entering the two generation processes. As a consequence α A a and α B a are not two simple rates but rather two (symmetrical) binomial distributions with respect to a: with a = 0, ..., Z describing the Z + 1 possible numbers of loaded A type FPs. It has to be noted that if the number of A type FPs in a unit size vesicles is given, the number of B type FPs is determined as b = Z − a, where Z describes the total number of FPs in a unit size vesicle.

Budding matrices {
The parameter η > 1 describes the increased loading affinity of an A type FP to its corresponding vesicle compared to the affinity of a B type FP to such a vesicle. Symmetrically, η characterises the preferential loading of coatB-based vesicles with B type FPs. As an example the { A γ ai−1 ai } matrix entries are calculated in the following. The probabilities P A and P B of loading a single A or B type FP, respectively, to a coatA vesicle which buds off from an organelle of size i with a i FPs are given as: These elementary probabilities enter the recursion formula for A γ ai (m, n) describing the successive loading of a vesicle with at present m FPs comprising n ≤ m FPs of type A: A γ ai (m, n) = A γ ai (m − 1, n − 1) · P A (a i − n + 1, iZ − a i + n − m) + A γ ai (m − 1, n) · P B (a i − n, iZ − a i + n + 1 − m), for m = 1, ..., Z, n = 0, ..., (a i − a i−1 ) and A γ ai (0, n) = 0 ∀n = 0, A γ ai (0, 0) = 1. The budding matrices in Eq. (7) are finally determined by A γ ai−1 ai = A γ ai (Z, a i − a i−1 ).