Dynamic balance between vesicle transport and microtubule growth enables neurite outgrowth

Whole cell responses involve multiple subcellular processes (SCPs). To understand how balance between SCPs controls the dynamics of whole cell responses we studied neurite outgrowth in rat primary cortical neurons in culture. We used a combination of dynamical models and experiments to understand the conditions that permitted growth at a specified velocity and when aberrant growth could lead to the formation of dystrophic bulbs. We hypothesized that dystrophic bulb formation is due to quantitative imbalances between SCPs. Simulations predict redundancies between lower level sibling SCPs within each type of high level SCP. In contrast, higher level SCPs, such as vesicle transport and exocytosis or microtubule growth characteristic of each type need to be strictly coordinated with each other and imbalances result in stalling of neurite outgrowth. From these simulations, we predicted the effect of changing the activities of SCPs involved in vesicle exocytosis or microtubule growth could lead to formation of dystrophic bulbs. siRNA ablation experiments verified these predictions. We conclude that whole cell dynamics requires balance between the higher-level SCPs involved and imbalances can terminate whole cell responses such as neurite outgrowth.


Introduction
Neurite outgrowth is an early event that changes the state of neurons [1] and allows neurons to develop axons and dendritic trees that connect to other neurons and become electrically active. When there is nerve injury and the axons are severed, the process of regeneration which is similar to neurite outgrowth often fails, resulting in the formation of bulbs at the ends of regenerating axons. These bulbs are called dystrophic bulbs [2]. Although some of the cellular pathways involved in axonal regeneration, failure and formation of the dystrophic bulbs are known [3,4] the mechanisms by which these lead to failure are not well understood. In this study we have used an integrated computational and experimental approach to understand the origins of dystrophic bulbs and more generally the conditions under which whole cell responses such as neurite outgrowth can be maintained. The central question we asked in this study is whether quantitative imbalances between subcellular processes can lead to the formation of dystrophic bulbs.
Whole cell responses that involve, both morphological and physiological changes, are complex because they can engage many subcellular processes (SCPs) that are together responsible for whole cell functions. SCPs include biochemical pathways and cell biological processes involved the functions of different organelles. Coordinated changes in the activities of SCPs often lead to a change in cell state, such as moving to a more differentiated phenotype. Detailed understanding of how different SCPs collaborate to produce an integrated whole cell response requires experiments, along with mathematical models that capture the dynamics of individual SCPs as well as the interrelationship between the dynamics of multiple SCPs. Based on the biological detail that is represented by the SCPs they can be classified as Level-1, Level-2 and Level-3 SCPs. Higher level SCPs (i.e. SCPs with smaller numbers) represent a cell biological function that combines the functions of their lower level SCP children [5]. To understand whole cell responses, it is necessary to delineate the quantitative coordination between the dynamics of the SCPs. In the case of neurite outgrowth, three Level-1 SCPs, Production of Membrane Components, Vesicle Transport & Exocytosis and Microtubule Growth need to be coordinated. Membrane components such as membrane lipids and proteins are synthesized in the cell body and transported to the growth cone via vesicle transport. Vesicles that fuse with the growth cone plasma membrane deliver their membrane lipids to the neurite tip, causing the neurite shaft to grow. Microtubule constitute the structural element of the neurite shaft, and they need to grow as the shaft increase in length. Tsaneva-Atanasova et al [6] developed a model of the interaction between microtubules and membrane vesicles in regulating neurite outgrowth. The framework of this model serves as a useful starting point for more extensive models where we can consider how a group of distinct but related SCPs such as Vesicle Budding, Vesicle Transport and Vesicle Fusion along with a second group of SCPs that involves Microtubule Nucleation, Microtubule Growth and Conversion of Dynamic to Stable Microtubules interact to enable the growth of neurites. Such detailed models could serve as basis for understanding how a single network can regulate the relationship between kinetic parameters and changes in the levels or activities of molecular components within the different SCPs. To develop a detailed model of neurite outgrowth based on interactions between SCPs of different levels, we constructed a multi compartment ordinary differential equation model to extend and integrate the microtubule growth model developed by Margolin et al. [7] with the vesicle transport model by Heinrich and Rapoport [8] (Heinrich and Rapoport, 2005). We have also incorporated more recent experimental data to make the model contemporary.
Many mathematical models of neurite outgrowth have been developed [9,10] and these models have provided valuable insights into how complex biological processes can be modeled at various levels of description. Often these models abstract the details of the underlying mechanisms and consequently it is difficult to decipher how the balance between the various SCPs is achieved to generate whole cell responses. Some studies have focused on specific facets of the neurite outgrowth process such the role of signaling and regulatory motifs in determining how particular neurites are selected to become axons [11]. Such models are useful in understanding the regulation of complex cellular responses. In developing our dynamical model of neurite outgrowth, we have made a critical change from our previous approaches to dynamical modeling in developing this model. We had typically used a bottom-up approach where individual reactions are assembled into modules that have emergent behaviors [12,13]. In contrast, the approach we use here can be called top-down as we made a starting postulate that the three Level-1 SCPs involved in neurite outgrowth need to be appropriately coordinated to produce a steady velocity of neurite outgrowth. We start with experimental observations of neurite outgrowth velocities and use numerical simulations to determine how their Level-2 children and Level-3 grandchildren SCPs need to be coordinated to allow neurite outgrowth with the observed velocities and under kinetic parameter constellations that match additional experimental data. Our simulations show redundancy between Level-3 sibling SCP activities. Changes in function of one Level-3 sibling SCP can be compensated by alteration of function in another. However, such generalized reciprocal relationships do not exist between Level-2 SCPs. Change in activity in one Level-2 SCP cannot be simply compensated by activity adaptation of one of its sibling Level-2 SCPs but requires adaptations of multiple sibling SCPs. Consequently, loss of coordinated whole cell function is more likely. A direct consequence of this loss of coordination during neurite outgrowth is the formation of dystrophic bulbs. We used the findings from analyses of the computational models to develop predictions regarding effects of modulating higher SCPs on whole cell function. We tested the predictions experimentally and found them to be valid.

Model description: Integration of subcellular processes involved in vesicle production, transport and exocytosis and microtubule dynamics to model neurite outgrowth
Neurite outgrowth and axonal regeneration share common SCPs such as vesicle production & transport and microtubule growth in an orderly fashion leading to a growing shaft and a growth cone at the tip as shown schematically in Fig 1A left panel. Failure of growth and the subsequent formation of dystrophic bulb (Fig 1A right panel) can be hypothesized as failure of SCPs to function "properly" with respect to each other. To understand how this lack of coordination between SCPs occurs it is necessary to have a detailed description of the SCPs and their relationships as part of setting up the computational model we provide a detailed description of SCPs and their relationships. Among the major SCPs that are involved in neurite outgrowth (NOG) are the Level-1 SCPs Production of Membrane Components, Vesicle Transport and Exocytosis and Microtubule Growth (Fig 1B). Each of these SCPs is composed of multiple Level-2 children SCPs whose activities need to be coordinated to ensure NOG without complication.  [35] as well as disorganized microtubules [3]. (B) The three major SCP types, Production of Membrane Components, Vesicle Transport & Exocytosis and Microtubule Growth, describe essential sub-cellular functions involved in neurite outgrowth. We labeled them as Level-1 SCPs and identified their Level-2 children and Level-3 grandchildren SCPs. (C) To model membrane production at the Trans-Golgi network (TGN) and delivery to the growth cone plasma membrane (GC-PM) via vesicle transport, we extended a dynamical model that simulates vesicle transport between the endoplasmic reticulum and the Golgi [8]. In our version, new synthesized membrane is added to the Trans-Golgi network from where it is transported to the growth cone plasma membrane by vesicular transport. Vesicles bud from the TGN, move along the microtubules via active kinesin mediated transport and fuse with the growth cone plasma membrane, leading to neurite shaft growth. The vesicles pass through three intermediate compartments: the cell body cytoplasm (CBC), the neurite shaft cytoplasm (NSC) and the growth cone cytoplasm (GCC). In each compartment kinesin has a different affinity for microtubules, leading to varying fractions of vesicles that are actively transported via kinesin along the MT. Retrograde vesicles are generated via endocytosis at the growth cone plasma membrane and move along the microtubule towards TGN through dynein mediated active transport. Vesicle budding at the TGN or GC-PM is mediated by the interaction of recruitment factors and coat proteins, and vesicle fusion with the TGN or the growth cone membrane by the formation of SNARE complexes between vesicle(v)-SNAREs and target(t)-SNAREs, catalyzed by local tethering machineries. Motor proteins are bound to vesicles through motor protein receptors. A G /B G label vesicles that bud from the TGN with coat protein A/B, A PM / B PM vesicles that bud from the GC-PM with coat protein A/B. Since at steady state almost all anterograde vesicles bud with coat protein B from the TGN, and almost all retrograde vesicles bud with coat protein A from the GC-PM, we highlighted their transport routes in brighter colors. (D) To simulate MT growth, we consider two different MT pools, stable and dynamicMTs. After nucleation new MTs are added to the pool of dynamic MTs that is characterized by alternating phases of growth and catastrophic breakdown. The frequency and duration of these phases depend on the tubulin concentration (regulated by the SCP Tubulin Sequestration) and GTP hydrolysis rate (SCP GTP Hydrolysis). Consequently, the length distribution of the dynamic MT pool and the degradation rate of dynamic MTs depends on the activity of both SCPs (S2 Fig). Dynamic MTs are either degraded or converted into stable MTs that form the growing neurite scaffold.
Since our major focus was on vesicle transport and microtubule growth, we further categorized Level-2 SCP function by introducing their Level-3 children SCP. SCPs with the same parent SCPs are defined to be siblings. This limited scope enabled us to develop detailed, yet manageable computational models. It is important to note that at this stage our model does not consider local actin dynamics at the growth cone, since we particularly focus on SCPs that are involved in extension of the neurite shaft; i.e. membrane lipid and protein synthesis and delivery and microtubule scaffold growth.
The level-1 SCP Production of Membrane Components has two Level-2 children SCPs, Membrane Lipid Production and Membrane Protein Production, that are incorporated into our model via constant synthesis rates of the individual components at the Trans-Golgi Network (TGN) (Fig 1C; Reaction 1 in S1 Table) [14,15].
The Level-2 children SCPs of the Level-1 SCP Vesicle Transport & Exocytosis are Vesicle Budding at Trans-Golgi Network and Vesicle Endocytosis at Growth Cone Plasma Membrane (GC-PM), Anterograde and Retrograde Microtubule-based Vesicle Transport, Vesicle Fusion with Trans-Golgi Network and Vesicle Exocytosis at Growth Cone Plasma Membrane. Each of these six Level-2 SCPs has two Level-3 children SCPs. To simulate the dynamics of the vesicle transport SCPs we extend a dynamical model of bidirectional membrane lipid and protein transport between the endoplasmic reticulum and the cis-Golgi [8]. In our model the two organelles are the TGN and the GC-PM ( Fig 1C). The level-2 SCP Vesicle Budding at the TGN consists of the Level-2 SCPs Coat Protein Recruitment and Formation at the TGN and Vesicle Invagination and Scission at TGN. The former is mediated by Recruitment Factor 1 that resides at the TGN, the latter is mediated by a cell body specific budding rate (Reaction 2 in S1 Table). Membrane proteins at the TGN are recruited into the budding vesicle via association with the coat proteins (S2 Table). The Level-2 SCP Anterograde Microtubule-Based Vesicle Transport & Exocytosis that moves vesicles from the cell body cytoplasm into the neurite shaft cytoplasm and from there into the growth cone cytoplasm involves the Level-3 SCPs Kinesin Recruitment to Vesicle and Kinesin Mediated Vesicle Transport along Microtubule. Kinesin motor proteins are recruited to the vesicle via kinesin receptors that are part of the vesicle membrane. Here, we assume that kinesin is available under saturating conditions, so that the kinesin concentration per vesicle is determined by the number of kinesin receptors. Based on the number of kinesin bound kinesin receptors we calculate the likelihood that a vesicle is bound to the microtubule (S1 Table, Reactions 6/7). This likelihood also depends on the affinity of kinesin to the microtubule, incorporated into our model as the fraction of Microtubule (MT)-bound kinesin. The fraction of MT-bound kinesin decreases from the cell body cytoplasm to the neurite shaft cytoplasm and growth cone cytoplasm (S3 Table). A low fraction of MT-bound kinesin in the neurite shaft cytoplasm ensures that only a few vesicles (~9%) are moving forward, in agreement with experimental data [16]. The large pool of stationary vesicles in the neurite shaft cytoplasm could serve as a membrane reservoir that allows quick recruitment of additional membrane to the GC-PM on demand. In the growth cone cytoplasm almost all vesicles are unbound, allowing their fusion with the GC-PM.
The Level-2 SCP Vesicle exocytosis at GC-PM is mediated by the interplay of the Level-3 SCPs Vesicle Tethering and Vesicle Fusion at GC-PM. The former is realized in our model via a growth cone specific tethering rate, the latter via complex formation between vesicle(v)-SNARE V and target(t)-SNARE Y that resides at the GC-PM (S1 Table, Reaction 8). Retrograde vesicle transport from the GC-PM to the TGN is incorporated by a similar set of Level-2 and Level-3 SCPs that consist of a symmetric set of membrane proteins and rate constants. We assume a high fraction of dynein is attached to the MT in the neurite shaft cytoplasm, so that most of the retrograde vesicles are actively transported towards the TGN, without formation of a vesicle reservoir.
Vesicles that fuse with the GC-PM, increase its surface area. In our model, we set the GC-PM surface area to 50 μm 2 based on experimental data that describe the size of the growth cone P-domain as about 50 μm 2 [17,18]. We use the size of the P-domain, since this is the area of vesicle fusion with the growth cone [19]. For simplicity, we still refer to this area as the GC-PM. Any membrane surface area that is added to the GC-PM and would increase its size beyond 50 μm 2 is added to the neurite shaft (S1 Table, Reaction 9). Membrane proteins are not added to the neurite shaft, since our model is based on the assumption that the diffusion of the membrane proteins into the neurite shaft is prevented by intra-membranous diffusion barriers or cortical cytoskeleton proteins [20] that ensure a highly specialized growth cone plasma membrane. Neurite shaft length arises from its membrane surface area (S1 Table, Reaction 10). For simplicity, we do not consider increase of the neurite diameter at this stage. The neurite shaft cytoplasm lies within the neurite shaft and grows with the lengthening neurite shaft. Most of the anterograde vesicles in the neurite shaft cytoplasm are not actively transported along the microtubule but are stationary components of the membrane reservoir. Consequently, the growing neurite shaft cytoplasm acts as a sink for vesicles and vesicle membrane proteins.
Microtubule (MT) growth is mediated by the interplay of the Level-2 SCPs Dynamic Microtubule Growth and Stable Microtubule Growth. [21,22,23,24]. Dynamic MTs are characterized by periods of MT growth and catastrophic breakdown, the latter can either lead to MT disappearance or rescue, followed by a new growth period. The time frame of these periods is in seconds or minutes. Stable MTs do not show such periodic growth behavior and we assume that there is no degradation of stable MTs. Dynamic MTs are continuously generated at a specified nucleation rate (Level-3 SCP Dynamic MT Nucleation) (Reaction 11 in S4 Table, Fig 1D). The length of dynamic MTs (Reaction 12, S4 Table) as well as the rate of MT catastrophic breakdown depends on the free tubulin concentration that is regulated by the SCP Tubulin Sequestration (Reaction 13, S4 Table) [25,26,27] and on the hydrolysis rate of tubulin-bound GTP (SCP GTP Hydrolysis). The latter SCP determines the size of the GTP cap at the tip of the dynamic MT that protects it from catastrophic breakdown [28].

Conceptual frame work for modeling whole cell dynamics arising from interactions between SCPs
Our model uses whole cell responses in a top-down based manner that started with the literature curation of SCPs involved and their assignment to different levels of biological detail ( Fig  2, Step 1a), as described above. Like the organization of our Molecular Biology of the Cell Ontology [5], Level-1 and Level-2 SCPs describe more general, Level-3 SCPs more specific cell biological or biochemical functions. SCPs of the different levels are in parent-child relationships. SCPs are represented as coupled biochemical reactions that are transformed into coupled ordinary differential equations (ODEs) (Fig 2, Steps 1b/c). We use experimental data to define model constraints that the simulations must match (Fig 2, Step 1d). For systematic analyses we develop analytical solutions for the prediction of kinetic parameters of SCPs at steady state (Fig 2, Step 2a). Here, we use the term steady-state to describe a continuous whole cell response without variation of amplitude. We use the analytical solutions to predict steady-state kinetic parameters and values that match the predefined model constraints (Fig 2, Step 2b). Through systematic parameter variation, we generate multiple steady-state solutions that allow for the classification of dependencies between SCPs at different levels (Fig 2, Step 2c). To investigate the effect of SCP dysfunction on whole-cell responses, we generate numerical simulations based on steady-state parameter sets after arbitrary modification of kinetic parameters of a single SCP followed by experimental verification (Fig 2, Step 2d).

Simulations, analytical analyses and predictions
To characterize the coordinated activities of these Level-3 SCPs, we ran an initial set of simulations of the growth of a single dynamic MT using a stochastic growth model [7] (S1 Fig). Since we simulate NOG over hours and days, while dynamic MT growth behavior changes in seconds and minutes, we do not simulate individual dynamic MTs, but the whole population of dynamic MTs. For this population simulation, we developed two formulas based on the results of the model by Margolin et al. They describe the dependence of the average growth length of dynamic MTs and of the catastrophic breakdown rate that leads to complete disappearance of the dynamic MT on the effective tubulin concentration and the GTP hydrolysis rate (Reactions 12/13 in S4 Table; Table, Reaction 15). To calculate the length of the microtubule scaffold within the neurite, we sum the length of all dynamic and stable MTs and divide it by the number of MTs per neurite cross-section (Reaction 16/17 in S4 Table).

Experimental determination of varying neurite growth rate velocities
In primary cultures, neurons put out neurites at varying rates. To quantify the variability and determine the velocities with which neurites grow, we used a high content imaging system (IN Cell Analyzer) to track neurite outgrowth in primary rat cortical neurons.

Fig 2. Top-down based SCP modelling of whole cell dynamics.
Our pipeline for the generation of dynamical models that allow the analysis of relationships between sub-cellular processes (SCPs) during whole cell responses consists of two major steps. First, we build multicompartmental ODE models for the simulation of SCP activities that are based on literature curated data (Step 1a-d). In the second step we develop analytical predictions for kinetic parameters at steady state that allow a systematic analysis of SCP interactions (Step 2a-d). https://doi.org/10.1371/journal.pcbi.1006877.g002 The experimentally observed velocities of neurite outgrowth were used as a top-down constraint for our model. Neurons were plated at a low density under optimal growth conditions on a permissive substrate of Poly-L Lysine. We incubated them for 16h and then documented neurite outgrowth every 6h for an additional 48h. For each neuron, we quantified the length of its longest neurite (the one most likely to become an axon). Representative images are shown in Fig 3A and 3(b), and a complete set of images is shown in S2 Fig. At each time point the neurites significantly differed in their lengths (Fig 3B). We identified the longest neurite for each neuron in a field (S2 Fig) and calculated their length distribution at every time point ( Fig  3C), identified the length of selected quantiles and calculated the growth velocity for these quantiles (Fig 3D). The neurite growth velocities ranged from 0 to~20 μm/h. We used these velocities as a top-down whole-cell response constraint in modeling the activities of the Level- 1, -2 and -3 SCPs. Increase in the whole-cell response, i.e. neurite outgrowth velocity, depends on the coordinated increase of membrane production and delivery of membrane vesicles to the growth cone as well as microtubule growth. Both types of SCPs need to facilitate growth with the same velocity, since they are interdependent. Vesicle transport is microtubule dependent and microtubule growth is limited by the length of the neurite shaft. Similarly, the increase in membrane production needs to be coordinated with increased vesicle transport. More vesicles would need to bud from the TGN, to be actively transported through the neurite shaft cytoplasm and to fuse with the GC-PM. These interdependencies form the basis for our top-down constraint.

Postulates governing SCP dynamics
In addition to the experimentally determined NOG velocities, we used experimental data from the literature for additional constraints at the level of SCPs. These constraints arise from a set of postulates of permissible SCP functions that are described in (S5 Table). We postulate that about~10% of the anterograde vesicles in the neurite shaft cytoplasm are bound to the microtubule and actively moving along the microtubule [16]. The other 90% constitute a membrane pool that can be recruited on demand to adapt to short-term increases in outgrowth velocity. Similarly, we assume that the vesicles that reside in the growth cone [29]are the first membrane reservoir for the quick recruitment of additional membrane and also define-in the same way -that 10% of the anterograde vesicles in the growth cone cytoplasm move, i.e. fuse with the growth cone membrane. This assumption is based on the experimental observation that GC membrane precursor vesicles reside in the GC cytoplasm before fusion with the membrane with a half-life of �14min [30], which would suggest that around 5% of the vesicles fuse with the membrane. Similarly, we assume that 90% of the retrograde vesicles in the neurite shaft cytoplasm and the cell body cytoplasm are moving. The total growth cone membrane is completely internalized within 30-60 min [31] which suggest an endocytosis rate between 0.833 μm 2 /min-1.66 μm 2 /min for a growth cone with a surface area of 50 μm 2 [17,18]. Based on these observations and the assumption that some of the endocytosed vesicles might be part of transcytotic re-arrangements during growth cone steering [32] or are transported to other growth cones and not to the cell body [33], we set the rate of membrane back transport from the growth cone to the TGN to 0.5 μm 2 /min. The backward transported membrane can also be called cycling membrane, since it continuously cycles between the TGN and the GC-PM without being added to the neurite shaft. For the microtubules we postulate that dynamic MTs should fill up 20 μm of the growing neurite. This is the length of the minor processes [34] that are characterized by continuous retraction and extension periods [1]. We assume that such periods are committed by continuous MT growth and retraction, so that the MT scaffold in these processes should mainly consist of dynamic MTs.
Using these constraints, we ran simulations to identify an initial solution at one intermediate rate of neurite outgrowth. From these initial simulations, we developed an overall analytical approach that allows us to predict the relationships that produce neurite outgrowth at a specified velocity.
We identified a set of parameters that generated neurite outgrowth with a velocity of 10 μm/h (Fig 4 and S3 Fig). Under steady-state conditions neurite, surface area continuously increased at a rate of 0.5236 μm 2 /min in agreement with the amount of membrane that is needed to increase the length of the growing neurite ( Fig 4A). The length of the microtubule bundles also increase proportionally to maintain the shaft growth rate (Fig 4B). For this we selected an effective tubulin concentration of 9 μM to achieve a length increase of the microtubule bundles with a velocity of 10 μm/h. The vesicles in the neurite shaft cytoplasm increased in parallel with a rate of~0.57 vesicles/min assuming a membrane area of 0.05 μm 2 per vesicle ( Fig 4C). In considering the values of the ordinate for the different SCPs (Fig 4D) it became clear that there are multiple combinations of kinetic parameters and amounts of components  Table)  Dynamic balance between vesicle transport and MT growth enables neurite outgrowth that allow neurite outgrowth at the same velocity without violating the constraints we imposed on dynamics of the SCPs.
To enable a systematic analysis of the redundancies within and between SCPs, we generated an analytical solution (Fig 5, and S4 Fig) for the prediction of kinetic parameters at steady state that enable neurite outgrowth at a given velocity without violating constraints on SCP dynamics (S5 Table). In our solution, we consider that an increase in neurite outgrowth depends on varying sets of protein amounts and kinetic parameters that allow increased net anterograde transport of membrane vesicles and continuous synthesis of new proteins and membrane to compensate for the accumulation within the growing neurite shaft cytoplasm. In addition to the postulated constraints on SCP dynamics, we pre-define independent protein amounts and kinetic parameters (S4 Fig, right side) to calculate the dependent protein amounts and kinetic parameters. To achieve a given NOG velocity a certain amount of lipid membrane needs to be synthesized and transported from the TGN to the GC. In addition to this membrane that will be added to the neurite shaft three other membrane types can be classified based on their Dynamic balance between vesicle transport and MT growth enables neurite outgrowth transport destination, allowing us to distinguish four membrane types (Fig 5A). These membrane types should not be confused with the four-different vesicle sets that transport the membrane between the TGN and GC-PM. Type I membrane buds from the TGN and is added to the vesicle reservoir in the growing neurite shaft cytoplasm. Type II membrane is transported from the TGN to the GC-PM and added to the neurite shaft as anterograde vesicles. Type III membrane is transported anterogradely from the TGN to the GC-PM, endocytosed at GC-PM and added to the NSC as retrograde vesicles. Type IV membrane is transported from TGN to GC-PM and from there back to TGN, and therefore is labeled cycling membrane. Initial and final anterograde and retrograde rates that describe budding and fusion processes at the TGN and GC-PM can be calculated by summing up the involved membrane types. The fluxes of the individual membrane types can be calculated based on the anticipated NOG velocity and model constraints (Fig 5B, Eq.(1)). To allow the fusion of sufficient anterograde vesicles with the GC-PM, a certain number of SNARE complexes need to be formed between v-SNAREs V and t-SNAREs Y (Fig 5B, (Eq. (2))). Consideration of the tethering rate allows the calculation of the needed amount and concentration of v-SNARE Vs. The concentration allows the determination of the initial and final forward fluxes of SNARE V (Fig 5B, (Eq. (3)). Similarly, concentrations and fluxes can be calculated for the v-SNARE U that mediates vesicle fusion at the TGN (Fig 5B, Eqs. (4) and (5) 6) and (7)). The production rates of v-SNARES (and any other membrane vesicle membrane proteins) can be calculated based on their vesicle concentrations and the number of vesicles that accumulate in the growing neurite shaft cytoplasm (Fig 5B, Eq. (8)).
MT growth at a given velocity (Fig 5C) depends on the conversion rate of dynamic into stable MTs (Fig 5D, Eq. (9)). By pre-simulation results using the model of Margolin et al, we obtained the dependence of the average length of dynamic MTs on the free (i.e. not sequestered) tubulin concentration and GTP hydrolysis rate. The average length allows the calculation of the number of dynamic MTs (Fig 5D, Eq. (10)), followed by calculation of the stabilization rate (Fig 5D, Eq.(11)) and subsequently the nucleation rate (Fig 5D, Eq.(12)).
The analytical solution was validated by comparing the analytically predicted outcome with the results of the numerical simulation for example sets of parameter and protein amounts (Fig 5E and S5 Fig). We used the predicted parameter and protein amount sets as starting points in our model and ran the numerical simulation for 5000 minutes (~3.5 days). Every 500 min we compared the anticipated model constraints with the numerically obtained values to estimate the accuracy of our analytical prediction. The analytical solution matched the anticipated model constraints with high accuracy, except the fusion rate for retrograde vesicles with the TGN that for was too high for the velocities 0 and 2.5 μm/h.

Complementary compensatory relationships between Level-3 sibling SCPs can contribute to robustness of response
We used the analytical solution to systematically analyze the relationship between different SCPs (Fig 1B). For our studies we selected growth velocities that covered the whole range of experimentally observed velocities (Fig 6). We predicted that Level-3 sibling SCPs would have complementary activities. An increase in the activity of one Level-3 SCP can be compensated by a decrease in the activity of at least one of its siblings Level-3 SCPs. To acquire a specified NOG velocity without violation of the model constraints there are multiple combinations of the activities of the two Level-3 SCPs Coat Recruitment and Formation and Vesicle Invagination and Scission (Fig 6A), the two children SCPs of the Level-2 SCP Vesicle budding at the The more a Level-2 SCP activity is generated by a Level-3 SCP that contains a vesicle membrane protein, the higher the SCP activity of membrane protein production to compensate for the loss of membrane protein in the growing NSC reservoir. (E) In dependence of the selected fraction of MT-bound kinesin, vesicles need a different amount of kinesin receptors to ensure NOG without violation of the model constraints. A higher kinesin receptor concentration per vesicle directly translates into the need for a higher TGN. Such SCP redundancies are also observed between the Level-3 SCPs Kinesin Recruitment to the Vesicle and Kinesin Mediated Vesicle Transport along the MT (Fig 6B) that define the activity of the Level-2 SCP Anterograde Microtubule-Based Vesicle Transport. Similarly, the SCPs Vesicle Tethering and Vesicle Fusion (Fig 6C), sub processes of the SCP Vesicle Exocytosis are characterized by an interaction that is based on an inverse relationship between their activities.
The balance of the Level-3 SCP activities to fulfill the Level-2 SCP function is not isolated at Level -3 but influences the activity of other SCPs as well. The more the activity of a Level-2 SCP arises from a Level-3 SCP that involves vesicle proteins, the higher the need for continuous membrane protein supply to compensate for the consumption of membrane proteins in the neurite shaft cytoplasm sink. This accounts for the kinesin receptor ( Fig 6E) as well as for the v-SNARE V (Fig 6F). Here, one might consider that the complementary SCP Kinesin-Mediated Vesicle Transport along the MT is determined by components along the whole microtubule while the complementary SCP Vesicle Tethering at GC-PM is determined by components of the tethering machinery that are locally restricted to the growth cone. The maintenance of the activity of the former SCP might therefore depend on a continuous synthesis of cytoplasmic protein that is not considered in our model. Similar compensatory relationships exist between Level-3 sibling SCPs involved in Microtubule Growth. The combined length of dynamic MTs that are available for conversion into stable MTs is regulated by an inverse relation between the activities of the Level-3 SCPs Dynamic MT Nucleation, GTP Hydrolysis and Tubulin Sequestration. The average dynamic MT length (Fig 7B) as well as the degradation rate of dynamic MTs (Fig 7C) depends upon the Level-3 SCPs Tubulin Sequestration (which regulates effective tubulin concentration) and GTP Hydrolysis. The balanced activities of these two Level-3 SCPs and the Level-3 SCP Dynamic Microtubule Nucleation (Fig 7D & 7E) determines the activity of the Level-2 SCP Dynamic Microtubule Growth. Similarly, to the vesicle transport Level-3 SCPs, any loss of function in one of these MT Level-3 SCPs can be compensated by a gain of function in one of the other.

Sibling Level 2 SCPS do not readily compensate for each other leading to dystrophic bulbs
The relationships between Level-2 SCP activities are more complicated. A decrease in the activity of one Level-2 SCP cannot simply be compensated by increase in the activity of one of its Level-2 sibling SCPs. To simulate mismatching activity of a Level-2 SCP on NOG, we first predicted a set of kinetic parameters that allow growth under physiological conditions as defined by our model constraints. Then we selected one of its Level-3 SCP children, followed by the perturbation of one kinetic parameter that belonged to this Level-3 SCP, so that this parameter no longer matched anymore to the predicted parameter set. The resulting unbalanced activities between the redundant Level-3 SCP siblings propagated to impaired Level-2 SCP activity that did not match with the activities of the other Level-2 SCPs. We then searched for adaptations in other parameters that allow restoration of neurite outgrowth under physiological conditions as defined by our model constraints. As shown for the example Level-2 SCP Vesicle Budding at the Trans-Golgi Network, multiple Level-2 SCP activities needed to be altered in a coordinated manner (Fig 8), e.g. Membrane Lipid and Protein Production kinesin production rate. Lines refer to pre-defined fractions of MT-bound kinesin: 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 (from top to bottom). (F) Different tethering rates demand different v-SNARE V concentrations per vesicle and consequently different V production rates to ensure NOG without violation of the model constraints (t-SNARE Y is kept constant). Lines refer to different tethering rates: 2.5 x10 -6 , 3 x10 -6 , 3.5 x10 -6 , 4 x10 -6 , 4.5 x10 -6 , 5 x10 -6 (from top to bottom).
https://doi.org/10.1371/journal.pcbi.1006877.g006 (simulated via reduced production rates), Anterograde Vesicle Transport along the Microtubule (simulated via reduced amount of Kinesin receptor) and Vesicle Exocytosis at the GC-PM (simulated via reduced amount of v-SNARE V). These activity changes were associated with a reduction of NOG velocity (from 10 μm/h to 5.5 μm/h), since the cycling rate was kept constant at 0.5 μm 2 /min.
To simulate the effect of mismatching Level-2 SCP activities on NOG, we reduced the activity of one Level-2 SCP (as described above, but without adaptation of any other Level-2 SCP activities). We selected one Vesicle Transport & Exocytosis SCP and one Microtubule Growth SCP for further investigation. To simulate the effect of the mismatching activity of the Level-2 SCP Vesicle Exocytosis at GC-PM, we selected the v-SNARE V that is involved in the activity of its Level-2 child SCP Vesicle Fusion at the Growth Cone. We predicted that a reduction in its activity does not influence NOG velocity, but overloads the growth cone with anterograde vesicles, a feature associated with dystrophic bulbs [35] (Fig 9A & 9B). To test our predictions, we selected the v-SNARE VAMP7 that is involved in vesicle fusion at the GC-PM [36]. Cortical neurons were plated into one side of microfluidic chambers. These chambers contain two sides that are separated by a microgroove wall. Small slits in the wall allow growing neurites but not cell bodies to cross into the other side. Consequently, the neuronal cell bodies are isolated from the growing axons. Two hours after plating the neurons we transfected the neuronal cell bodies with siRNA against VAMP7 or MTCL1. To allow the siRNA to have an effect on protein expression we incubated the transfected neurons for 48h before we performed an axotomy to eliminate all basal axon growth on the axonal side of the chambers. We waited another 48hrs to allow axonal regrowth under knock down conditions. Neurite outgrowth was quantified after immunostaining against β-III tubulin. Knock down success was documented by RT-PCR against VAMP7 or MTCL1(Suppl. Fig 6). In agreement with our predictions, VAMP7 knock down (KD) did not significantly decrease NOG velocity, as shown by the total length of neurites 48h after axotonomy (Fig 9C & 9D). Its main effect was the generation of dystrophic bulbs (Fig 9C & 9E).
To simulate impaired activity of the Level-2 SCP Stable Microtubule Growth we selected its Level-3 SCP child Microtubule Stabilization and perturbed the kinetic parameter stabilization rate. Our simulations predicted a significant decrease in the total MTB length that affects NOG (Fig 9F). We validated our predictions by knocking down the microtubule crosslinking protein MTCL1 that mediates MT stabilization [37]as described above (see S6 Fig for documentation of knock down efficiency). MTCL1 KO significantly reduced NOG velocity and increased the number of dystrophic bulbs (Fig 9G, 9H & 9I).

Whole cell responses depend on subcellular activities at different levels
Neurite outgrowth depends on the coordinated activities of various sub-cellular processes (SCPs). Of the different SCPs involved in our model, the types of Level-1 SCPs are critical. Three major sub-cellular functions that are covered by these SCPs are Production of Membrane Components, Vesicle Transport & Exocytosis and Microtubule Growth. We curated SCPs that are involved in these Level-1 sub-cellular functions and classified them as Level-2 and Level-3 SCPs that are in parent-child relationships. Based on the hierarchical relationships between these SCPs, we developed a dynamical model that incorporates the three types of Level-1 SCPs and documented that the activity of the SCPs involved need to be highly coordinated to allow neurite outgrowth at a certain velocity. To develop this whole cell model of neurite outgrowth, we merged two models [7,8] that contain biochemically well characterized and estimated parameters. We updated these models by inclusion of further reactions. However currently there is no suitable experimental framework to reliably estimate hundreds of unknown parameters for different reactions in mesoscale dynamical models such as ours. Hence our combined neurite outgrowth model does have several parameters that are unknown. Some of these were obtained from the published biochemical literature and are specified in Table S7. The remaining parameters are constrained by top down analytical solution approach we use. To make the model more manageable, our simulations were done with the following assumptions: (1) the size of the TGN & GC-PM (50 μm 2 ) is independent of the outgrowth velocity; (2) the rate of back transported membrane from the growth cone to the TGN is 0.5 μm 2 /min; (3) 90% of the anterograde vesicles in the neurite shaft cytoplasm and the growth cone cytoplasm serve as membrane reservoir (i.e. these vesicles are not actively transported along the MT or fusing with the GC-PM) and (4) 20 μm of the neurite microtubule scaffold consists of dynamic MTs and the rest of stable MTs. All these assumptions are based on previous experimental data and allowed us to identify the key features of the balance between SCP types required for a dynamic whole cell response. We developed an analytical prediction for steady state dynamics at the

Complementary dependence in lower level SCPs enables robustness of response, while imbalances in higher level SCPs are not readily compensated
We could show that the activities of Level-3 SCP siblings needed to generate the function of their Level-2 SCP parent inversely depends on each other. This applied to both Level-3 SCPs involved in Vesicular Transport & Exocytosis as well as Level-3 SCPs involved in Microtubule Growth. This complementary compensatory capability will confer to the cell the ability to mount a response in a robust manner.
In contrast, the activities of Level-2 SCPs must be highly coordinated to maintain the Level-1 SCP function. The impaired activity of one Level-2 SCP cannot be simply compensated by an increased activity of another Level-2 SCP. Dependencies between Level-2 SCPs are more complex. In case of Vesicle Transport & Exocytosis, loss of function of one Level-2 SCP might only be compensated, if multiple other Level-2 SCPs also reduce their function. Impaired Level-2 SCP activities were simulated by reduction of the activity of one of its Level-3 children SCPs without compensatory activity modification of the other.
Our predictions proposed that reduced exocytosis does not influence NOG velocity. We predicted that this results in an excess of anterograde vesicle in the GC that can be one of the reasons for the formation of dystrophic bulbs. In agreement with our prediction, knockdown of v-SNARE VAMP7 did not change NOG velocity but rather, it creates a pathological phenotype with multiple dystrophic bulbs instead of healthy growth cones.
A reduction in the growth of stable MTs was predicted to cause a reduction in NOG length, and it was experimentally confirmed by knock down of MTCL1. As the activities of the SCPs related to Vesicle Transport & Exocytosis are unchanged, in this case the reduction in neurite outgrowth will be accompanied by the appearance of dystrophic bulbs. This prediction was also experimentally confirmed.

Different subcellular processes can give rise to similar whole cell phenotypes
The observations in this study both computational and experimental demonstrate the critical need for balance and coordination between different types of subcellular processes required to mount a whole cell response. It is noteworthy that robustness of whole cell responses arises from the relative redundancy of lower level subcellular processes within each type. However at least for neurite outgrowth, such redundancies do not occur at higher level SCPs where, when imbalance occurs the whole cell function shuts down even though each high level SCP might be operational for a certain time period. It is the uncoordinated functioning of the higher level VAMP7. P1 rat cortical neurons transfected with siRNA against VAMP7 or scrambled siRNAs were stained for β-III tubulin to document neurite growth across the microgroove of the chamber. In agreement with our prediction, VAMP7 knock down did not decrease total outgrowth length. It caused a significant number of physiological growth cones to turn into dystrophic bulbs. (D) Analysis of quantified neurite lengths confirmed that there was no significant difference in outgrowth length. Shown are average values and standard deviations (three independent experiments). (E) Neurites growing from neurons after VAMP7 knock down showed a significant increase in the number of dystrophic bulbs, indicative for impaired NOG. (F) To simulate the effect of reduced MT stabilization, we reduced the stabilization rate that turns dynamic into stable MTs to 25%, 50% and 75% of the original prediction, without changing the predicted values for all other kinetic parameters. Numerical simulations propose that decreasing stabilization leads to a significant decrease in NOG outgrowth over the investigated time. (G) We validated our predictions by knock down of the microtubule crosslinking protein MTCL1. MTCL1 knock down caused reduced NOG outgrowth and dystrophic bulb formation. (H) Outgrowth quantification confirms significant decrease in NOG (3 independent experiments). (I) Mtcl1 also significantly increased the number of dystrophic bulbs. � p<0.05, �� p<0.01, ��� p<0.001based on unpaired t-test for all panels. https://doi.org/10.1371/journal.pcbi.1006877.g009 Dynamic balance between vesicle transport and MT growth enables neurite outgrowth SCPs that result in the morphological abnormalities like the dystrophic bulbs. Since lack of coordination can arise due to dysfunction of different higher level SCPs, we get similar dystrophic bulbs when we partially ablate a microtubule stabilizing protein or a vesicle fusion protein. This provides a striking example of phenotypic convergence in spite of molecular divergence.
At this stage our model only includes microtubule dynamics and vesicle transport. The inclusion of additional subcellular functions into the model, such as actin filament dynamics at the growth cone, could further document multiple molecular states that underlie the same phenotype. Actin filaments interact with microtubules at the growth and could be an additional factor regulating dynamic microtubule growth apart from the free tubulin concentration and GTP hydrolysis rate that our model considers so far. The same might account for the interaction between actin filaments and vesicles that plays an important role within vesicle endocytosis and exocytosis.
Overall the major conclusions of this study are that 1) elements of both robustness and fragility coexist within subcellular processes needed to produce whole cell responses; 2) since multiple types of subcellular processes are involved in mounting a whole cell responses changes in in different molecular entities or subcellular processes can result in the same aberrant phenotype; 3) quantitative imbalances between different types of subcellular processes are sufficient for alteration in cellular dynamics causing abnormal phenotype.

Modeling
Our dynamical model of neurite outgrowth is based on the integration of vesicle transport model of Heinrich and Rappaport [8] and microtubule dynamics model of Margolin et al [7]. Since both models are several years we updated the models with additional details as needed. The model contains a large number of variables and it is not possible to estimate all of the parameters from the experimental literature. Hence,we adopted top down modelling approach that is based a top level experimental constraint of velocity of neurite outgrowth that we measured and multiple molecular kinetic parameters that we obtained from previous models and the experimental literature. At the top-level, we made a dynamic steady state assumption: that to maintain a given growth velocity the rate of microtubule growth and membrane vesicle dynamics should be balanced. With this assumption, we estimated th needed parameter values for different experimentally measured neurite outgrowth velocities using an analytical solution approach which also helped us to predict relationships between sub cellular processes at steady state. A full description of the computational model is in S1 Text. The reactions are in S1 Table. The parameters are S2, S3, S4, S5 and S7 Tables and the descriptions of symbols are in S6 Table   Experiments Rat cortical neurons were prepared as described in Supporting Information Text. All procedures complied with IUCAC rules of Mount Sinai School of Medicine.
To measure the velocities of neurite outgrowth primary neuronal cultures were plated on 96well plates incubated for 16hrs. Then brightfield images acquired every six hrs. up to 70hrs using the IN Cell Analyzer. Neurite lengths were quantified using the Metamorph Image Analysis Software.
For the siRNA ablation experiments rat cortical neurons transfected with control and siRNA mixture of interest were grown in microfludic chambers and after axotomy neurite outgrowth was measured by confocal microscopy on a LSM880 microscope and quantified using Image J software. Effectiveness of siRNA based ablation was measured by changes in cognate mRNA levels by quantitative RT PCR. For details see S1 Text Supporting information S1 Fig. Dynamic model of microtubule growth. (A) Using the model by Margolin et al [35], we generated four growth profiles of single dynamic microtubules over a time period of 26,000 secs under varying effective tubulin concentrations (ranging from 8 to 10 μM) as well as varying GTP hydrolysis rates (ranging from 0.65/sec to 0.75/sec). Four example profiles are shown for (effective tubulin = 9 μM and hydrolysis rate = 0.67 sec -1 , 0.72 sec -1 , figures shown only initial 6000 seconds). (B) We identified the length distribution for each growth. Shown are the distribution of the four examples in (A). (C) The average length of the dynamic microtubule calculated for given range of effective tubulin concentrations, hydrolysis rates and used third order cubic polynomials to fit the data. We used all average lengths to generate a formula that describes the dependence of the average MT length on the effective tubulin concentration and the GTP hydrolysis rate. This formula was incorporated into the main model (D) To determine the degradation rate of dynamic MTs in dependence of the effective tubulin concentration and GTP hydrolysis rate, we counted how often the simulated dynamic MT of each growth profile falls below the length of 4 tubulin dimers that we defined as the threshold for complete catastrophic breakdown. We used the results to generate a formula that describes the dependence of the degradation rates on the effective tubulin concentration and the GTP hydrolysis rate. This formula was incorporated into the main model. (TIF)

S2 Fig. Neurites grow with different outgrowth velocities.
Neurons were dissected from rat cortical brain, plated on 96 well plates and incubated for 16h to allow initial growth, followed by Image acquisition every 6h up to 70h after plating. Images were pseudo colored and subjected to neurite length quantification via Metamorph.  Fig 4C that shows the number of vesicles). B G refers to anterograde vesicles that bud from the TGN with the coat protein B, A G refers to anterograde vesicles that bud from the TGN with the coat A, B PM refers to retrograde vesicles that bud from the GC-PM with coat protein B and A PM to retrograde vesicles that bud from the GC with coat protein A. Green lines refer to TGN or growth cone plasma membrane (GC-PM), orange lines refer to anterograde moving vesicles (B G , A G ) and blue lines to retrograde moving vesicles (B PM , A PM ). Solid lines represent those vesicles that are mainly responsible for membrane transport in the indicated direction (i.e. B G in the anterograde direction and A PM in the retrograde direction).  Table. Dissociation constants of membrane proteins for interactions with the coat proteins. Dissociation constants were selected to ensure preferred binding of proteins that are involved in anterograde transport to coat B and proteins that are involved in retrograde transport to coat A. High dissociation constants for stationary proteins for those coat proteins that would incorporate them into vesicles involved in back transport were selected to ensure that these proteins stay at their anticipated organelle. Orange: Proteins involved in anterograde movement, Blue: Proteins involved in retrograde movement, Italics: Proteins that cycle between the TGN and GC, Standard fonts: stationary molecules. k d values are taken from [19] which is normalized. (DOCX) S3 Table. Percentages of microtubule (MT) bound/unbound motor protein in the different compartments. A high fraction of MT-bound kinesin in the cell body cytoplasm ensures that vesicles are immediately transported forward after budding from the TGN. A low percentage of bound kinesin in the neurite shaft cytoplasm causes about 90% of the anterograde vesicles in the neurite shaft cytoplasm to be stationary and constitute a membrane mobilization reservoir that can be recruited upon demand. This fraction depends on the total amount of kinesin receptors (Fig 6B). No kinesin is bound to the MT in the GC cytoplasm, so that all anterograde vesicles are available for fusion with the GC membrane. The percentages of MT-bound dynein in the GC and CB cytoplasm were selected for the same reason. Since we do not assume that there is a membrane reservoir for retrograde vesicle, we selected a high percentage of bound dynein in the neurite shaft cytoplasm. (DOCX) S4  Table. Parameters and initial values. The parameters set allows growth at 10 μm/h as shown in (Fig 4).