A Computational Analysis Framework for Molecular Cell Dynamics: Case-Study of Exocytosis

One difficulty in conducting biologically meaningful dynamic analysis at the systems biology level is that in vivo system regulation is complex. Meanwhile, many kinetic rates are unknown, making global system analysis intractable in practice. In this article, we demonstrate a computational pipeline to help solve this problem, using the exocytotic process as an example. Exocytosis is an essential process in all eukaryotic cells that allows communication in cells through vesicles that contain a wide range of intracellular molecules. During this process a set of proteins called SNAREs acts as an engine in this vesicle-membrane fusion, by forming four-helical bundle complex between (membrane) target-specific and vesicle-specific SNAREs. As expected, the regulatory network for exocytosis is very complex. Based on the current understanding of the protein-protein interaction network related to exocytosis, we mathematically formulated the whole system, by the ordinary differential equations (ODE). We then applied a mathematical approach (called inverse problem) to estimating the kinetic parameters in the fundamental subsystem (without regulation) from limited in vitro experimental data, which fit well with the reports by the conventional assay. These estimates allowed us to conduct an efficient stability analysis under a specified parameter space for the exocytotic process with or without regulation. Finally, we discuss the potential of this approach to explain experimental observations and to make testable hypotheses for further experimentation.

Although experimental studies have provided invaluable insights for the underlying exocytosis mechanisms, the process of exocytosis is a typical example to show the difficulty in conducting an analysis at the systems biology level [24][25][26][27][28]. That is, while the biochemical reaction chain is straightforward and simple, the regulation in vivo of the system is complex. As many kinetic rates are unknown, and concentrations of proteins, complexes and substrates keep changing in both vivo and vitro environments, a biologically meaningful, global system analysis is intractable in practice. Earlier, Mezer et al. [10] proposed a computational platform to model the exocytotic process. They formulated these protein interactions into a sequential (feed-forward only without any regulation) interaction pathway to describe the exocytotic system dynamics.
In this paper, we utilize the exocytotic process as a model system to present a computational framework for system modeling and analysis. Similar to [10], we model the dynamics and architecture of the complex system by the ordinary differential equations (ODEs). First, we model the whole system by taking the regulatory elements into account. Second, we use a math techniques called inverse problem to estimate the rate parameters for the basic steps of biochemical reactions. Through the method, we are able to recover and optimize these parameters based on limited in vitro experimental data. Third, based on the above estimates, we can therefore approximately study the stability behavior of this system with and without MUNC18 regulation. We then attempt to explain experimental observations about different fusion efficiency caused by the change of SNARE proteins' concentration and multiple complexes in the SNARE-induced membrane fusion. Moreover, we make a few interesting predictions that can be verified by further experimentations.

Results and Discussion
The Protein Interaction Network of Exocytosis From the view of gene network, the exocytotic process is a sophisticated combination of sequential interactions of welldefined proteins and protein complexes [1][2][3][4][5][6][7][8][9][10][11][12][13]. As shown in Fig.1, it has three major components. The first step of the basic reaction component includes two membrane proteins, SNAP25 (synaptosome-associated protein, 25 kDa) and syntaxin, together forming the so-called tSNARE; here t means target, the plasma membrane where the vesicle is heading for. Another important protein is vesical-associated membrane protein (VAMP2), belonging to the category of vSNARE (vesicle). In the second step, the protein complex formed by tSNARE and vSNARE is the fundamental step for the membrane fusion. In our study we consider two regulatory components, which are MUNC18mediated and Ca 2z -dependent regulation pathways, respectively. On the other hand, from the view of systems biology the mechanism of this exocytotic process is a dynamics system capturing the temporal change of the concentrations of proteins and intermediate complexes, which can be formulated based on an ODE dynamic system, as shown below in details.
The basic steps. The well-known foundations [4,3,8] for this exocytotic processes are the following two reactions.
where the protein complex FHC stands for the four-helical bundle. Formation of FHC complex is the main step to promote membrane fusion, an essential part of exocytosis. In addition, there are several follow-up complex modifications. For instance, the function of complexin is as a clamp [20], resulting in where FC is the generic notation for the protein complex of FHC and complexin. Ca 2z -dependent regulation. Ca 2z is the main trigger for the initiation of intracellular exocytosis [4,14,15,5]. Suggested by [16,4], the regulation of Ca 2z is executed through stimulating synaptotagmin. A well-known mechanism is that Ca 2z binds with the SNARE complexes (FHC) and stimulates the fusion [2]. The reaction equations to characterize the mechanism regarding Ca 2z and synaptotagmin are given by where CaS stands for the complex of synaptotagmin and one Ca 2z ion, and where the generic notation Tsc represents the protein complex of tSNARE and synaptotagmin binding four Ca 2z ions, and FHC Ã represents the complex of FHC and synaptotagmin binding four Ca 2z ions. MUNC18-dependent regulation. MUNC18 is an important regulatory protein for the exocytotic system [6,17,22], through two different modes: (i) MUCNC18 associates with syntaxin to remove them from the assembly into the SNARE complexe at the beginning stage; and (ii) MUNC18 stimulates the fusion process by associating with FHC. These two reaction mechanisms can be written as follows. Figure 1. The whole process of fusion used in the mathematical model is shown. One direction arrows and symbol of 00 z 00 represent the reaction between proteins, ions and complexes, while full direction arrows connect two parts of a single reaction. Modified from [1][2][3][4][5][6]. doi:10.1371/journal.pone.0038699.g001 where FHC ÃÃ is the complex of MUNC18 and FHC to help the fusion process; and where Smc is the generic name for the protein complex of MUNC18 and syntaxin.

Mathematical Modeling for the Whole Exocytotic System
Putting together, we have formulated a mathematical model by the ordinary differential equations (ODE) to capture how the concentrations of different proteins and complexes vary with time and how they interact each other. Based on the law of mass action and Michaelis-Menten Kinetics, and using the conventional notation ½ : for the concentration, the ODE system is given by.
One may raise the question, due to the complexity of this network, whether we have empirical evidence enough to show the concept we try to put forward. In a recent article, we [29] have conducted a comparative network motif analysis for the Sec1/ Munc18-SNARE regulatory mechanisms through a comprehensive compile of experimental data from different species and different cell types. In spite of some differences in details that have been shown important for cell-specific and species-specific system behaviors, we confidently conclude that Eq.(1) may conceptually represent the basic dynamic system that is likely universal. Some comments about Eq.(1) are presented below.
Dynamics of Ca 2z . The dynamics for Ca 2z in exocytosis is complex. To be analytically feasible, we assume that during the fusion process, concentration of Ca 2z ions at active zone is temporal dependent. Thus, the dynamics of Ca 2z ions around the region of fusion (active zone) can be characterized as.
where S(Ca 2z ) is the recruitment source of calcium. Nevertheless, the in vivo concentration of Ca 2z ions may stay at a roughly constant level as both external and internal sources may have kept the balance of Ca 2z . In this case, Eq.(2) can be replaced by the simplest form Ca 2z~C onstant. Self-association of syntaxin. We notice that self-association of syntaxin is possible such that i ' Ei=E-i syntaxin i , where i~5,6,7,8, syntaxin i represents the complexes made of i syntaxins [12]. Hence, if we take the effect of self-association into account, the system of Eq.(1) needs to be modified as follows: we have the equation for the concentration of syntaxin.
and additional four equations to describe the dynamics of selfassociated complexes, that is, where i~5,6,7,8 Mass conservation. One can show that the ODE system of Eq.(1) complies with the detailed balance principle and the mass conservation. For instance, because the only products of the reactions involving MUNC18 are Smc and FHC ÃÃ , the change of concentration of MUNC18 is only relevant to the concentrations of these two complexes. Indeed, for the subsystem that only involves MUNC18, Smc and FHC ÃÃ , we obtain.
since there are no Smc and FHC ÃÃ initially. Spatial effect. Denote all of the variables (concentrations) in Eqs.(1)-(5) by a vector U so that the ODE system can be rewritten in a concise form of dU=dt~f (U), where f is vector of functions on the right hand side of each equation. If the spatial effects of proteins and protein complexes are considered, this system should be generally written as follows.
where D U stands for the vector of diffusion coefficients of proteins, complexes and ions. Study of reaction diffusion equations described by Eq.(6) would be interesting particularly for the problems related to the developmental process.

Estimation of Reaction Rate Parameters
The whole system for the exocytotic process as described in Eq.(1) is a typical example to show the general difficulty in systems biology [24][25][26][27][28]. While the biochemical reaction chain is simple, the regulation in vivo of the system can be very complex. In addition, most paremeters remain unknown in this ODE system, including the initial concentrations of different proteins and complexes, and the reaction rates in both vivo and vitro environments. Hence, it is almost intractable in practice to carry out a global system analysis. As a first step to overcome this difficulty, we attempt to estimate the rate parameters for the basic steps of exocytotic process. Among different methods, we choose the technique of inverse problem that has two advantages: the required data size is small, and the algorithm guarantees the uniqueness and efficiency [25].
The fundamental subsystems. As the well-known machinery, chemical reactions.
are fundamental for membrane fusion. The behavior of this subsystem involving only proteins SNAP25, syntaxin and VAMP2 can be described as Using the inverse problem technique [30] to estimate rate parameters requires initial concentrations. In the following the symbol ½u(0) is used for the concentration of variable u at time t~0. From the experimental data [9], we set the initial condition for system Eq.(10) to be: ½SNAP25(0)~½Syntaxin(0)~9½VA MP2(0)~9mm=L, and ½tSNARE(0)~½FHC(0)~0. It should be noticed that the estimation of kinetic rate parameters are usually insensitive to the initial conditions, as verified by our simulation studies (not shown).
Reparameterization for data-fitting. In the experimentation, researches use the fluorescence intensity, x(t) to measure the time-dependent fusion process. The relationship between the concentrations of core complexes (FHC, FHC Ã , FHC ÃÃ ) and the fluorescence intensity, x(t)~C(½FHC,½FHC Ã ,½FHC ÃÃ ), needs to be addressed in some details. Some experimental studies such as [11] suggested that the function C can be roughly considered to be linear when the signal strength is far below the saturated level. In this case we have.
Estimation by the inverse problem algorithm. To recover the appropriate reaction rates, we apply technique introduced by [27] to Eq. (10). Some useful theorems are presented in the section of Materials and Methods. Using the data from [9], the identified parameters are shown in the table 1. We compare the numerical results based on the identified parameters with experimental data in Fig.2, and the error is e 10 {4 .

Stability Analysis of the Fundamental Subsystem
Estimation of rate parameters of the subsystem Eq.(10), as summarized in Table 1, allows us to carry out the stabilizing analysis under a specified parameter space. Considering the subsystem Eq.(10) with v 2 instead of x(t), we first study the fundamental subsystem without any regulation, under the initial concentrations u 1 (0), u 2 (0), u 3 (0), v 1 (0) and v 2 (0) for proteins SNAP25, Syntaxin, and VAMP2, and protein complexes tSNARE and FHC, respectively. While the formal mathematical treatment is shown in the section of Data and Methods, below we discuss about the biological interpretations.
Our analysis has shown that the final steady state level of the fusion is highly dependent on initial concentrations. Obviously, three proteins (SNAP25, syntaxin and VAMP2) must exist at t~0 so that u 1 (0)w0, u 2 (0)w0 and u 3 (0)w0. It is reasonable to assume no any fusion (here measured by FHC) at the initial time point, which means v 2 (0)~0. The only case we have to deal with carefully is the initial concentration of tSNARE complex, v 2 (0). This is because in vivo, tSNARE is already preformed in the plasmic membrane; and then carried by vesicles, vSNARE (VAMP2 in our case) binds with it to generate fusion. In this sense, we assume v 1 (0) §0 in general.
(A ) Case-A assumes that the initial concentration of SNAP25 and syntaxin are the same such that c 1~0 . Denote K~k {1 =k 1 , provided k {2 vvk 2 , we have shown there are two locally stable-steady states, denoted by P 1 and P 2 , respectively, corresponding to c 2 wc 3 or c 2 vc 3 . If c 2~c3 , the degenerated steady state P is also stable. Since we are mostly interested in the final steady state-level of fusion, i.e., v v 2~½ FHC, the biological meaning of above stability analyses can be summarized in Table 2. In short, for the system involving SNAP25, Syntaxin and VAMP2, we should only consider two types of initial conditions: If the initial conditions only include initial concentrations of SNAP25, Syntaxin and VAMP2, but no tSNARE, the final steady state of fusion ½FHC is equal to the least initial concentration of SNAP25, Syntaxin and VAMP2. In the case of non-zero initial concentration of tSNARE (½tSNARE 0 w0), however, the final steady state ½FHC can be much higher as long as the initial concentration of VAMP2 (vSNARE) is sufficiently large. This case is particularly interested because in vivo, SNAP25 and Syntaxins may have been already preincubation (preformed) into tSNAREs on the plasmic membrane, before vSNARE proteins (VAMP2 in our case) approach, as carried by vesicles.
Our analysis explains why the outcome of fusion process depends on the way to put these three proteins into the system [21]. One is the sequential process: SNAP25, Syntaxin and VAMP2 proteins are added into the system in order such that virtually no tSNARE protein complex has been formed when the reaction begins. The other one is the preformed process: After SNAP25 and Syntaxin proteins have been preincubation (preformed) into tSNARE, VAMP2 proteins are then added to initiate the fusion reaction. Numerical simulations have shown that the preformed process reaches the steady state much faster than the sequential one (Fig.3), which is consistent with in vitro experimental data (the embedded panel) [21].

Stability Analysis on MUNC18-dependent Regulation
We furthermore study the stability behavior of the system involving the regulatory protein MUNC18. As discussed above, the MUNC18-dependent regulation has two types: (i) It binds tightly to a closed conformation of sytanxin that precludes the syntaxin's involvement in the fusion process, suggesting that MUNC18 inhibits fusion by regulating the formation of tSNARE. And (ii) it can assemble with SNARE complexes (FHC) to accelerate membrane fusion in late stages when the concentration of four helical bundles (FHC) is high enough.
We make the following assumptions to simplify the subsystem with MUNC18-dependent regulation. Considering the situation that tSNARE has been preformed and reaction of SNAP25, sytanxin and tSNARE has reached the equilibrium, we claim that the function of MUNC18 can be characterized as follows.
where FHC ÃÃ is the complex of MUNC18 and FHC, and it behaves similar to FHC to help the fusion process; and tSNAREzMUNC18 ? p fSNARE{MUNC18g where p stands for the binding rate of MUNC18 onto the syntaxin in closed conformation. In the above reactions, the concentrations of four helical bundles, FHC, and four helical bundles binding with MUNC18, FHC** reflect the level of fusion. It has been shown that the disassociation rate of MUNC18-syntaxin complex is very small comparing to the binding rate, so that the second reaction is considered as an irreversible one. Introducing variables ½vSNARE~u 3 , ½tSNARE~v 1 , ½FHC~v 2 , ½MUNC18~u 4 , and ½FHC Ã Ã~v 3 , the subsystem involving MUNC18 is rewritten as.
Using the mathematical approaches similar to the case of no regulation, we have studied the stability of system Eq. (11). Assume the binding rate p, reaction rates k +4 are in the range given by the reference [17][18][19], we have shown the existence of steady states of Eq.(11), including bi-stability. As the result has been rigorously presented in the section of Data and Methods, we are mainly interested in the final fusion level, as measured by F~FHCzFHC ÃÃ . Under the assumption that the initial concentration of FHC** is zero, i.e., v 3 (0)~0, we interpret our results as follows.
(i) If ½MUNC18 0 {½FHC 0 ƒ½tSNARE 0 ƒ½VAMP2 0 z ½MUNC18 0 , there exist two bi-stable states for the final fusion levels: One is the high fusion level, which is given by.
In this case, at the steady state, the concentrations of free MUNC18 and free vSNARE (VAMP2) are virtually zero, which mean all of these proteins exist in the form of FHC and/or FHC**. The second steady-state is the low fusion level (F low ), the   up-bound of F low is actually F z low~F high , whereas the low-bound is given by.
On the other hand, the low steady state fusion level, F low is somewhere between F { low and F high . In this case, the steady state, the concentrations of free MUNC18 and free tSNARE are virtually zero, which mean all of these proteins exist in the forms of FHC and FHC**.
(ii) Otherwise, there exists only one steady state that is locally stable, and the final fusion F is somewhere between (0,F high ).
Hence, with the regulation of MUNC18, the steady states of final level of fusion is controlled by the initial concentration of MUNC18: The behavior of bi-stability exists only when the initial concentration of MUNC18 is intermediate, whereas the boundary is determined by initial concentrations of tSNARE, VAMP2 and FHC. A lower or higher ½MUNC18 0 results in a single steady state of the final fusion level. Moreover, the final fusion level depends on ½tSNARE 0 , suggesting that the preincubation (preform) of tSNARE is an important factor. Indeed, using numerical simulations, we have shown that for the system involving SNARE proteins, complexes and MUNC18, preformed assays have two advantages over the sequential one: first, preincubation advances reaction rates; second, preincubation support more fusion than sequential assays. Finally, we comment that the regulation mechanism of MUNC18 may be threshold dependent, i.e. there exists an optimal threshold ? which depends on the initial concentration of tSNARE and four helical bundle only, such that MUNC18's regulation function during fusion is maximized when the initial concentration of MUNC18 reaches the threshold (Xia et al, unpublished results).

Conclusive Remarks
In this study, we present a framework for modeling protein interaction network which are involved exocytotic process. The framework is based on classic chemical kinetic model that generates insights into system dynamics and stability. The computational experiments and mathematical analysis reveal that the frame reconstruct biological experimental observation successfully and is able to provide useful predictions.

Simulation Procedures
The kinetics simulation and analysis of the whole system or the subsystems were implemented in Matlab7.0R. Differential equations were solved using the ODE23s routine. For testing the robustness of parameters, we generated 2000 random parameter sets using Latin Hypercube Sampling when all parameters are varied +30% relative to their original values, with a a uniform distribution for each parameter.
The concentrations of reactant proteins are given in molar units. For non-soluble proteins such as vSNARE and VAMP2, we followed the work in [10] and based the protein concentration estimation on the concentration of secretory vesicles in molar. During the exocytotic process, the size of vesicle pools varies with respect to different cell types from 200 to 3000. Hence, the molar concentration of vesicles was estimated in the range of 0.2-30 nm. Accordingly, the concentration of VAMP2 is considered to be in an identical range of vesicle concentration (0.2-30 nm) [10]. The tSNARE proteins such as SNAP25 and syntaxin are thought to be vastly expressed in vivo and the studies [1][2][3][4][5][6] evaluated the concentration of these protein in a range of 0.1-100 mm. The essential regulatory protein Munc18 is known to be expressed at much lower levels, compared to SNARE proteins, with the concentrations in range of 1-30 nm [3][4]10].

Algorithm for the Estimation of Rate Parameters
To recover the appropriate reaction rates, we apply technique of solving the inverse problem introduced by [27] to Eq. (10). Some useful results are presented below. To be concise, the ODE system Eq.(10) is written as A(p)U~0, U 0 is the initial conditions, and the parameter set p~(k +1 ,k +2 ,c).
The inverse problem claims that the parameter identification of Eq.(10) is equivalent to the optimization problem of.
subject to A(p)U~0 and U(0)~U 0 , where P is the parameter space in R 5 z and J(p) is regularized energy functional where M and a are Tikhonov regularization parameter [25], Q is parameter-data mapping, U e is experimental data, and e DD{pDD 2 is the penalty function to guarantee the positivity of reaction rates. The rational of the inverse problem is based on the following theorem: Suppose the solution of Eq.(10) U~(u 1 ,u 2 ,u 3 ,v 1 ,x) is smooth, where ½0,T) is the observation time. Then, given observed data on each time point in ½0,T), the parameters identified by the inverse problem are locally unique with respect to the initial condition. Under the assumption that all of the reaction rates are roughly constant, the optimization problem is solved through a gradient-based method. The brief algorithm is sketched below: 1. Given initial condition U 0 , solving ODE system (3:1) by fourth order RK and mapping it on the observation data set, 2. Gradient representation: using forward difference to approximate +J, 3. Applying steepest descent to approach the global minimum starting with some initial guess, 4. Using adjoint scheme to approximate Hessian + 2 J of J, 5. Using the approximate solution given by step 2 as initial guess , and using Quasi-Newton method with + 2 J to find the appropriate parameter.
Case-A. Assume the initial concentrations of SNAP25 and syntaxin are the same so that c 1~0 . Denote K~k {1 =k 1 , provided k {2 vvk 2 . There are two stable steady states, P 1 and P 2 : (i) If c 2 wc 3 , we have.
Case if c 1 §0 a n d c 3 §c 2 zc 1 , t h e s t e a d y s t a t e i s P 3~( c 1 ,0,0,c 3 {c 2 {c 1 ,c 2 ) and it is stable locally; (iv ) if c 1 §0, c 3 ƒc 2 zc 1 , and c 3 §c 1 , the steady state is P 4~( c 1 ,0,c 2 zc 1 {c 3 ,0,c 3 {c 1 ) and it is a stable node locally.
Case-C. A more general case is that the reaction ratios where Note that P 1 ,P 2 and P 3 are locally stable nodes. When K is small sufficiently comparing to the concentrations of SNARE proteins and complexes, P 3 is reduced to.
Proof. A concise proof is presented below. From the definition of c 1 to c 3 , straightforward calculation simplifies the fundamental subsystem Eq.(10) as follows.
For Case-A that c 1~0 , Eq. (22) can be further simplified to be.
Notice that reaction rates recovered from the experimental data imply k {2 =k 2 e 10 {10 M, so that compared to the concentrations of SNARE complexes, k {2 =k 2 is negligible. Thus, the v 2 -nullcline determined by Eq.(23) so that. Denote K~k {1 =k 1 , u 1 -nullcline is given by.
The steady states are yielded by intersecting the nullclines, and the biological interesting steady states are therefore given by P 1 and P 2 . Straightforward calculation implies the steady states P 1 and P 2 are a pair of opposite vertexes, and the relationship of c 3 and c 2 determines the choice of these two steady states. If c 2 wc 3 , the only possible steady state is P 1 ; if c 2 vc 3 , the only possible steady state is P 2 .
To investigate the stability of those steady states, we calculated the corresponding Jacobian for system Eq.(23) and then evaluate the two eigenvalues, denoted by l 1 and l 2 , respectively. For steady state P 2 , two eigenvalues for the Jacobian have no zero real part because of c 2 vc 3 , so that steady state P 2 is a hyperbolic point of system Eq. (23). By Hartman-Grobman theorem, there exists a homeomorphism mapping the trajectories of Eq. (23) in an open set containing P 2 onto trajectories of its linearized system in an open set containing P 2 . Furthermore, the homeomorphism preserves the parameterizations by time. Therefore, local behaviors of P 2 is characterized by its corresponding Jacobian, leading to l 1 v0 and l 2 v2. Therefore, steady state P 2 is stable locally.
For steady state P 1 , we calculated the corresponding Jacobian and showed that none of the eigenvalues has zero real part, so that P 1 is a hyperbolic steady point of system Eq. (23). The local behavior of trajectories of (Eq.(23) in the neighborhood of P 1 is characterized by its linearized system with respect to P 1 . As the trace of correspondiong Jacobian is less than zero, we show P 1 is stable locally.
In the same manner, we have shown the results presented in case-B and case-C.
Proof: The proof is similar to the case of system Eq.(10) without regulation.