Nonequilibrium description of de novo biogenesis and transport through Golgi-like cisternae

A central issue in cell biology is the physico-chemical basis of organelle biogenesis in intracellular trafficking pathways, its most impressive manifestation being the biogenesis of Golgi cisternae. At a basic level, such morphologically and chemically distinct compartments should arise from an interplay between the molecular transport and chemical maturation. Here, we formulate analytically tractable, minimalist models, that incorporate this interplay between transport and chemical progression in physical space, and explore the conditions for de novo biogenesis of distinct cisternae. We propose new quantitative measures that can discriminate between the various models of transport in a qualitative manner–this includes measures of the dynamics in steady state and the dynamical response to perturbations of the kind amenable to live-cell imaging.

ii. Chemical conversion of an A particle to a B particle at site i with rate uf [m A i ] that depends on the "mass" (number of particles) m A i of A through the flux-kernel f [m A i ]. Similarly, a B →C conversion can occur at site i with rate vf [m B i ], and so on 1 . iii. Fission,movement and fusion of single particles with species-dependent rates: An A particle fissions from site i, moves in the anterograde (or retrograde) direction with rate w A f m A i (or w A f m A i ), and fuses with the mass on the neighboring site. Similarly a B (or C) particle fissions and moves forward with rates w B f m B i (or w C f m C i ) and backward with rates w B f m B i (or w C f m C i ). iv. Breakage of a finite fraction (α) of an aggregate and movement in the anterograde direction with rate D (α=1 corresponds to cisternal progression, and α<1 to sub-cisternal movement. ).
v. Exit of full aggregates or single particles from boundaries, with the same rates as those for movement in the bulk.
The vectorial nature of transport is modeled through the asymmetry in the anterograde/retrograde particle movement rates, and is parametrized as: Fig. 1 of main paper). The asymmetry factors γ A , γ B , γ C for the three species can be different, representing differential degrees of recycling of A, B, C, . . . particles to the ER (a generalization of [1]).
The parameter space of this model is quite large, encompassing the injection rate a, the interconversion rates u, v, . . . the fission rates w A , w B , w C , . . . and the corresponding asymmetry factors γ A , γ B , γ C , . . ., the aggregate movement rate D, the breakage fraction α, and finally the form of the function f itself.

S1.2 Rationalizing the model
Below we comment in some detail on various aspects of the model and also highlight its strengths and weaknesses as a model of Golgi biogenesis. 1. Representing the three-dimensional Golgi by a one-dimensional model: Our model is a spatial model that explicitly incorporates distance from the cis end as a relevant variable.
In the interest of analytical tractability, we include only one spatial dimension corresponding to the cis to trans direction, which is the main direction of molecular traffic, and also the axis along which the Golgi exhibits biochemical polarity. The 1D model can thus be considered an effective model obtained by integrating over the two directions perpendicular to this cis to trans direction. However, a 1D model of this sort cannot address questions related to the shape of individual cisternae or if trafficking pathways are branched. 2. A,B,C particles represent molecules in different stages of processing in the Golgi: In constructing this model, we imagine an A particle to be the equivalent of the unprocessed protein molecules arriving at the cis-Golgi from the ER, and particles of type B, C, D as representing proteins in different (successive) stages of processing. Likewise, an aggregate which is primarily of type B is analogous to a cisterna with a large fraction of semi-processed molecules in an early stage of processing.
Thus, in our model, the processing stages 'A', 'B', 'C' of the primary constituent molecules of a cisterna define the chemical identity of the cisterna. 3. Modeling enzyme-mediated cargo modifications as simple Poisson processes: The A→B and B →C conversions which we have treated as Poisson rate processes in the model are catalyzed by enzymes in the real system. A more realistic model of the Golgi could include both A, B, C 'cargo species' and E A , E B 'enzyme species', with specific enzymes having an affinity for specific cargo types. This kind of a detailed enzyme plus cargo model would allow us to model feedbacks wherein the distribution of B particles influences the distribution of the corresponding E B enzymes, and is in turn influenced by E B molecules which activate the B →C conversion. Thus, in general, we expect self-organization in this extended model to be more complex. Nevertheless, the effective rate-based model studied in this paper incorporates the sequential interconversion process using just a few parameters, and has the advantage that it is more amenable to analysis than a complex model. Moreover, as a spatial model that incorporates both chemical and transport processes in the Golgi, it provides a valuable framework for constructing a more detailed model with enzymes. 4. Model parameters are composites of many biophysical rates: The various particle exchange and interconversion rates in our model are effective or composite rates. For example, the particle movement rate combines both the rate at which a vesicle buds from an aggregate, as well as the rate at which it fuses with the next aggregate. Similarly, the B →C modification rate is determined jointly by the rates of attachment and detachment of enzyme molecules with B molecules, the concentration of the enzyme molecules, the reaction rate between the enzyme and the B molecule etc. . Because of the composite nature of the rates, it is not straightforward to find the corresponding experimentally measurable parameter. However, a model of this sort with a relatively small number of composite parameters allows for an economical description of the system, which makes it easier to identify and elucidate the qualitative effects at play in the system.

S1.3 Dynamical equations for mass
The following dynamical equations describe the time evolution of the average mass (number of particles) of each species at each site i in the lattice, in accordance with the elementary moves allowed in the model: where . . . indicates averaging over ensembles. The mass of each species at each site changes due to single particle exchange and aggregate movement between neighboring sites, as well as interconversion at any given site. If the system attains steady state, then the time derivatives in eqs. (S1a)-(S1c) can be set to zero, so that the equations take the form of flux-balance conditions for each site.
S2 Analysis of the pure Vesicular Transport (VT) case S2.1 Insights from the analytical solution In the pure VT limit (with D=0), eq. (S1) can be solved in the steady state, by setting time derivatives to zero and taking a continuum limit i/L→x in space, which transforms the equations into a set of coupled, second-order ODEs. Solving these ODEs yields For the three-species model, the solutions are quite involved, and it is more convenient to solve eq. (S1) numerically. Figure S1 Figure S1: (a) Equation (S1) can be solved numerically in the continuum limit to obtain =0.66, u=0.01125, v=0.00194, msat =200, Ksat =14.14). (b) The approximate m A (x) , m B (x) , m C (x) profiles (dashed lines) derived using eq. (S2b) are reasonably close to the actual mass profiles (solid lines) from simulations.
The key characteristics of the mass profiles can be illustrated using a two-species version of the model with 'A' and 'B' particles that undergo A→B conversion, and move with the same directional bias γ A =γ B =γ. In this case, the analytical steady state solution for f [m A (x)] and f [m B (x)] is relatively simple [1]: Below, we discuss some salient features of this solution:  Parameter combinations in the blue region lead to runaway growth of B particles; parameters in red region correspond to steady state, and the bold line is the boundary between the two phases. For parameters in the shaded red region, m B >2000 at the peak of the B profile (as deduced from the approximation in eq. (S2b)).
parameter combinations in the blue and red regions lead to runaway growth of B particles and steady mass of B particles respectively. Physically, runaway growth arises because the outgoing vesicular fluxes at a site saturates if the mass at the site is much higher than m sat . Thus, at sufficiently high rates of influx, the outgoing flux from a location may fail to balance the incoming flux, leading to runaway growth. A more detailed analysis shows that typically, runaway growth occurs only in specific regions of the system, even as the rest of the system attains steady state [1].

2.
Obtaining sharp peaks with large but finite mass: The transition from steady state to runaway growth at a particular site is associated with divergence of the average mass at that site (see eq. (S2b), where the mean mass m diverges as f [m] →K sat ). Thus, by tuning the system to be close to the phase boundary between steady state and runaway growth (represented in fig. S2 by the solid black line), the peak concentration of particles for any species can be made very high. (For instance, see fig. S2, where parameter combinations in the shaded region result in a steady state in which the peak concentration of B is greater than 2000 particles at a site). Moreover, even though the f [m] profile is relatively smooth, the m profile becomes very sharply peaked when the system is poised close to the phase boundary, leading to a sharply localized B-rich region. 3. Tuning the locations of mass peaks: The location of the peak of any mass profile (here the m B (x) profile) is determined by two length scales, l c and l t , which govern the rise and fall of the f [m B (x)] and m B (x) profiles near the cis and trans ends respectively. For the two-species case, eq. (S3) can be analyzed (see [1]) to show that (i) the length scale l c is small if the ratio u/w A is large: For high values of u/w A , most A particles undergo an A→B conversion before they can travel into the bulk, so that the concentration of B particles rises steeply close to the source itself, leading to a small l c (ii) l c is large if the asymmetry factor γ is high: As γ increases, the proportion of B particles traveling in the anterograde direction rises, shifting the region of high B concentration in this direction, thus leading to large l c . (iii) The length scale l t , which governs the spatial variation of f [m B (x)] near the trans end, only depends on the asymmetry factor γ, and is large when γ is small: The presence of the particle 'sink' at the trans end lowers particle concentration in this region. If there is extensive recycling or movement of particles in the retrograde direction (corresponding to γ close to 1/2), then the effect of the sink is 'transmitted' in this direction, lowering the local mass farther away from the trans end, resulting in large l t . The locations of the maxima of mass profiles in the 3-species model are also determined in a similar manner by the ratio of interconversion to fission rates and the degree of anterograde-retrograde asymmetry of particle movement. 4. Coarsegraining: From eq. (S3), we can see that the mass profile remains unchanged under the scaling: L→λL, u→λ 2 u, (γ −1/2)→λ −1 (γ −1/2). This provides a basis for coarsegraining the system, allowing us to reduce the model to a small number of 'effective' lattice sites. However, for the coarsegraining procedure to be meaningful, this effective number should still be large enough that the continuum approximation used in deriving (S3) holds.

S2.2 Structure formation with different types of flux kernels f [m]
The above analysis suggests that the flux kernel f [m Z ], which encapsulates how the vesicle fission and conversion rates depend on the mass of the parent aggregate, is a crucial determinant of structure formation.

I. Mass-independent rates:
III. Michelis-Menten type of rates: The above forms of f [m Z ] emerge naturally for processes that are catalyzed by enzymes: If the number of enzyme molecules is much lower than the number of A, B, C particles available for fission or conversion, then the reaction rate is enzyme-limited and roughly independent of m A , m B . . . [case (I)]. If the enzyme concentration is in excess of the reactant particles, then the reaction is substrate-limited and the reaction rate is proportional to the number of A, B, C particles available for reaction, leading to  fig. S2]. As discussed in sec. S2.1, we can solve for f [m Z (x)] in steady state, and approximately derive the mass profiles m Z (x) from these solutions, using the following approximations: m Z rates emerge if we assume that the mass at any location in the 1D model is obtained by integrating over the 2D disc-like structure. A 2D structure with m particles typically has √ m molecules at its perimeter. Assuming that only the perimeter molecules are available for reaction, we get reaction rates that are proportional to Equation (S5a) has been derived explicitly for the single-species model without interconversion in [2] and appears to hold in the presence of interconversion as well in the limit f   Width of maxima: In order to obtain non-overlapping compartments, i.e., well-separated peaks in the total mass profile, the width of the maxima must be much smaller than the distance between adjacent maxima. Sharp peaks (with small width) can be obtained only for flux kernels f [m] of types (I) and (III), that saturate at large m. As evident from eqs. (S5a) and (S5c) S4.

Sensitivity of mass profiles to small variations in rates:
Changes in model parameters can alter the peak locations of m Z (x) profiles as well as peak concentrations. However, in the regime of small interconversion rates (u, v w A , w B ) that we consider here, the locations are relatively insensitive to changes in parameters, and perturbations primarily alter the peak concentrations of some or all m Z profiles.
Below we consider how mass profiles respond to a change in influx, by measuring [δm Z / m Z ]/[δa/a] which is the ratio of the relative change δm Z / m Z that occurs at a location with average mass m Z , to the corresponding change δa/a in influx.  Figure S6: Structure formation with pure VT and MM kinetics for fission and chemical conversion but with θ=1/2 for fission and θ=1 for conversion. Structures are qualitatively similar to the case with the same value of θ for both processes (figs. 4(e), 4(f)). Parameters: a=1, D=0, wA = 0.124691, wB =0.049383, wC =0.040741, γA =0.5, γB =0.64, γC =0.724, u=0.001062, v=0.000136 This analysis also sheds light on the degree of fine tuning required to generate compartments. Since the mass profiles are most (least) sensitive to small changes in parameters when f [m] is of type I (type II), it follows that the maximum (minimum) degree of fine tuning is required in this case. By the same token, an intermediate degree of fine tuning is required for structure formation with type III f [m].
So far we have considered identical flux control functions for both fission and chemical conversion. However, it is quite plausible that the two processes are governed by different kinetics, for instance if both obey MM kinetics but with different saturation scales m sat and exponents θ in eq. (S2a). Figure S6 depicts structure formation in the pure VT limit, allowing for this possibility, by assuming θ=1/2 for fission (as would be the case if vesicles break off only from the rims of cisternae) and θ=1 for chemical conversion (assuming that Golgi-resident enzymes are distributed throughout the cisterna). While the analytical solutions described above are no longer valid for this case, structure formation and even dynamics is qualitatively similar to the less general scenario (with identical θ and m sat ) considered in the main paper.

S2.3 Pure VT model with higher number of species
The pure VT model can be generalized to include more species of particles; here we consider a 4-species version with A→B →C →D sequential interconversion. Figure S7(a) shows the average mass profiles for a set of parameters chosen such that the maxima of the B , C and D profiles are sufficiently well separated. The maximum of each mass profile can be made sharp by tuning the maximum of the corresponding f [m Z ] profile to be close to the saturation value K sat , so that the total mass profile has three distinct peaks  5, γB =0.63, γC =0.66, u=0.00875, v=0.00146, msat =200, Ksat =14.14, m sat =100, K sat =10.
corresponding to morphologically separate B, C and D compartments. The model can be extended in a similar fashion to include more number of species, and thus, obtain a higher number of compartments.

S2.4 Effect of homotypic fusion
We use numerical simulations to study the effect of introducing homotypic fusion in the pure VT model. Homotypic fusion refers to the tendency of an A (or B) vesicle to fuse with a neighboring aggregate with a rate that increases with the number of A (or B) vesicles in the target aggregate. To implement homotypic fusion in our model, we assume that an A particle at site i to hops to its (forward) neighbor j with a rate given by: The first bracketed term represents enzyme-mediated fission from site i with MM kinetics. The second term represents homotypic fusion of the fissioned A particle with the mass at site j at a rate that increases with the mass of A particles present at the site, but saturates to a constant value beyond a certain mass m sat . The hopping rates for B particles can be written in a similar fashion. We consider a scenario where in addition to vesicular transport, a finite fraction α of an aggregate is allowed to break off and move forward with rate D. This is akin to the Cisternal Progenitor model of [5], and is referred to the VT-dominated limit in the main paper. Breakage and movement of cisternal fragments at a non-zero rate D is sufficient to ensure that there is no runaway growth, even for very large influx rates a.
To illustrate this point we analyze a very simple case where the system consists of a single site (L=1), and where there is no interconversion. The average mass on the site evolves as: For a/w A >K sat and D=0, the incoming and outgoing flux at the site necessarily fail to balance, as f [m A ] cannot exceed K sat , leading to runaway growth of mass. Now consider the case with non-zero D. As m A increases, the average outgoing flux (which includes the term αD m A ) also increases, until it balances the incoming flux a. At this point, d m A /dt must become zero, ensuring that m A is constant at long times. Thus, as long as the outgoing flux has a component that keeps increasing with the mass m at the site, there can be no runaway growth.
In order to preserve the well-separated compartments that emerge in the pure VT model, the rate D must be smaller than the single particle fission rates. As D increases, the peaks in the mass profile become shallower, until adjacent peaks can no longer be resolved ( fig. S8). Thus, we consider VT-dominated scenarios where cisternal movement acts as a secondary track to vesicular transport. This eliminates runaway growth while maintaining sharp peaks in the mass profile (see also figs. 2(c) and 2(c ) of the main paper).

S4 Typical steady state configurations in aggregate representation
Steady state configurations of mass for various transport models (see fig. 2 of the main paper) can also be shown in an alternative 'aggregate representation' where the local mass of A, B, C at spatial position x is represented by the heights of the green, red, blue columns at x. This representation is particularly useful in tracking changes in mass due to perturbations, for instance in movies S9-S12, or for detecting fluctuations of the mass profiles about the time-averaged concentrations (see solid black lines in fig. S9, also fig. S4).

S5 Dynamical measurements
We consider below two kinds of dynamical measurement which address: (i) how the number of tagged particles decays in time after the full system is tagged at t=0, and (ii) how the compositional entropy changes during reassembly of the system from an unpolarized state.

S5.1 Dynamics of (fluorescently) tagged particles
We perform measurements inspired by iFRAP (inverse Fluorescence Recovery After Photobleaching) experiments in which the cargo pool in the entire Golgi region is fluorescently tagged and highlighted and then the subsequent decay of fluorescence monitored to infer the exit kinetics of the cargo molecules [3]. The observation of exponential decay of fluorescence in these experiments has been used to argue against the cisternal progression hypothesis [3]. To further explicate these findings, we implement a similar protocol in our model, by tagging all the particles present in the system at an arbitrary time instant t=0 in steady state, and monitoring how the proportion of tagged particles decays in time. Figure S10(a) shows M tag (t)/M tag (0) (averaged over 100 ensembles) for each of the four cases -pure VT, pure CP, VT-dominated and CP-dominated. A clear linear decay in tagged particle concentration occurs only for the pure CP case. Exponential decay is observed for both the pure VT case and the VT-dominated case where cisternal movement provides a secondary track for traffic. Further, even in the CP-dominated case, where the bulk of the transport is through cisternal progression, there is a significant deviation from linear decay. This suggests that the observation of exponential decay in fluorescence experiments is consistent with cisternal progression as long as there is concomitant vesicular movement. A similar conclusion is reported in [4] where the iFRAP experiments in [3] were quantitatively analyzed to show that the exponential decay can be explained by a number of alternative scenarios involving vesicular efflux in addition to cisternal movement.

S5.2 Dynamics of reconstitution after disassembly
We consider below the dynamics of re-assembly of the system in the pure VT and pure CP scenarios, and in the process also distinguish reassembly from de novo formation. In our terminology, reassembly refers to the reconstitution of the system from an initial state in which the molecules and markers that comprise  the cisternae are all present in the Golgi region but distributed in a uniform, unpolarized manner. De novo formation, on the other hand, refers to the re-formation of structures starting from an initially 'empty' state, i.e. with molecules resorbed into the ER and/or dispersed far into the cytoplasm.
To study re-assembly in Monte Carlo simulations, an initial state is constructed by redistributing all the A,B,C particles (that constitute the steady state structures for a specific set of rates) uniformly across the system. We then switch on all transport and interconversion processes at the same rates and monitor how the number of A,B,C particles in different regions evolves in time (see also movies (S7) and (S8)). As in the case of de novo formation, this information can be used to compute the region-wise compositional entropy S Ω (t) (see main paper for a formal definition). As the structures re-assemble, S Ω (t) decreases with time in the pure VT case and changes non-monotonically in the pure CP cases for the region 0. 5<x≤0.75 (fig. S10(b)). In general, S Ω (t) can show complex non-monotonic behaviour in both limits depending on the region of interest. This is in contrast with de novo formation where the dynamics of the region-wise entropy provides a clear way of discriminating between the two mechanisms, increasing with t during formation in the pure CP case, and decreasing with t in the pure VT case (see fig. 4(f) of the main paper).
The contrast in the dynamics of S Ω (t) during de novo formation and reassembly in the pure CP limit highlights the sensitivity of dynamical measurements to the initial state of the system. The difference between the two dynamics can be rationalized as follows. During de novo formation, early compartments mature into medial and then trans compartments, thus creating mixed compartments, and increasing the compositional entropy. Eventually as this process goes on, the creation of mixed compartments due to maturation of early compartments is balanced by the loss of mixed compartments due to maturation into late compartments, so that S Ω reaches a time-independent value. However, in the early stages of de novo formation, when this balance has not yet been achieved, S Ω (t) increases with time. In the case of reassembly, on the other hand, early as well as late components are already present in the system in a completely random, unpolarized manner. Reassembly thus essentially involves the re-emergence of polarity in the system, and is thus qualitatively different from de novo formation.
shift significantly from its unperturbed state (solid black line). The proportionate change is higher near the peaks of the mass profiles than at the valleys. Different compartments respond over appreciably different time scales, with the cis cisternae attaining their new size much earlier than the trans cisternae. The systematic change in cisternal size is easily distinguishable from steady state fluctuations. (S11)-(S12): A 20% drop (S11) or rise (S12) in influx causes a modest change in cisternal size. Dashed and solid black lines represent the typical cisternal size (obtained by averaging over many ensembles) in the unperturbed state and after a 20% drop/rise in influx respectively. Note how the cisternae in the movies can be significantly larger or smaller than this average size making it difficult to distinguish the systematic change in size due to altered influx from usual stochastic fluctuations.
Response to exit block [(S13)-(S14)] Movies (S13)-(S14) show how cisternae respond to a blockade of the exit site at the trans end in the pure VT (S13) and pure CP (S14) limits. In both limits, exit blocks induce a piling up of particles at the exit site, with little or no change in the rest of the system. In our model, this pile-up continues indefinitely in time, thus pointing towards the need for an additional mechanism that stabilizes the enlarged trans cisterna in the presence of an exit block.