Kinetic stability analysis of protein assembly on the center manifold around the critical point.

BACKGROUND: Non-linear kinetic analysis is a useful method for illustration of the dynamic behavior of cellular biological systems. To date, center manifold theory (CMT) has not been sufficiently applied for stability analysis of biological systems. The aim of this study is to demonstrate the application of CMT to kinetic analysis of protein assembly and disassembly, and to propose a novel framework for nonlinear multi-parametric analysis. We propose a protein assembly model with nonlinear kinetics provided by the fluctuation in monomer concentrations during their diffusion.
RESULTS: When the diffusion process of a monomer is self-limited to give kinetics non-linearity, numerical simulations suggest the probability that the assembly and disassembly oscillate near the critical point. We applied CMT to kinetic analysis of the center manifold around the critical point in detail, and successfully demonstrated bifurcation around the critical point, which explained the observed oscillation.
CONCLUSIONS: The stability kinetics of the present model based on CMT illustrates a unique feature of protein assembly, namely non-linear behavior. Our findings are expected to provide methodology for analysis of biological systems.

Background

Numerical simulation based upon multi-parametric kinetic equations is the principal methodology for the analysis of the behavior of biological systems. Researchers often encounter a number of parameters in the governing equations of the system. Here, we introduce the center manifold theory (CMT) for simplification of the study of dynamic biological systems. CMT provides mathematical prescription for carrying out reduction of the number of parameters near the steady state, as well as information regarding the stability of the steady state. As a result, simulation is oriented to illustrate behavior around the critical point, at which system behavior drastically changes in the qualitative structure. The observable change is termed bifurcation, and the threshold values of the parameters are referred to as critical values or bifurcation values. The aim of this study was to provide a simple algorithm for the application of CMT to multi-parametric kinetic equations, in order to clearly illustrate the behavior of the biological system. The CMT has been applied to the Lotka-Volterra model of predator–prey system to provide important simulation results [1, 2]. In addition, several pioneering studies have applied CMT to neural network analysis [3]. Time-delay and diffusive effects play important roles in bifurcation phenomena [1, 4]. However, to date, there are few applications of the CMT to biochemical reaction models. We previously reported a model of cell signaling systems using non-linear kinetics and demonstrated the phase transition phenomenon via a numerical simulation [5].

Pivotal protein-protein interactions during cytoskeleton formation were selected as the application model for the present CMT method. Among the interactions between protein monomers, tubulin and actin polymerization are well-known events that have been analyzed using the numerical method [6-10]. The physical robustness of the cytoskeleton is based on the biophysical properties of actin and tubulin. In particular, various mathematical models have been proposed to explain the kinetic behavior of tubulin assembly [6-11]. A theory of polymerization of macromolecules has been established on the basis of the kinetic model of aggregation [12, 13]. Oosawa and Asakura previously reported that polymerization is similar to micelle formation or crystallization, and that there is a critical monomer concentration above which monomers effectively polymerize. The authors additionally suggested that the nucleation step represents the rate-limiting step for polymerization. Nucleation and growth occur in parallel during the progression of polymerization. There is a gap in free energy change between initial nucleation and progression of linear polymerization [13]. The stable nucleus for polymerization consists of trimers or tetramers, and the growth of aggregates through elongation/ dissociation follows the formation of a thermodynamically unfavorable size of the nucleus. In the current study, we focused on the polymerization in the absence of de novo nucleation and the interaction between polymer and monomer (PM) interaction.

For stable growth, the lifespan of tubules is controlled by a guanidine triphosphate (GTP)-cap that forms at their ends [14]. The structure and motility of growing tubules is influenced by intrapolymeric Brownian motion and fluctuation; this provides elasticity to the microtubules [15]. Polymerization/de-polymerization is controlled by binding of adenosine triphosphate (ATP)/GTP, resulting in the assembly of monomeric proteins. The intermittent transition between slow growth and rapid shrinkage in polymeric assemblies of microtubules is termed dynamic instability [14]. Numerous models have been proposed to explain this instability; in particular, Zapperi and Mahadevan successfully identified two parameters: a structural mechanical parameter that characterizes the ratio of longitudinal to lateral interactions in an assembly, and a kinetic parameter that characterizes the ratio of timescales for growth and conformation change. These parameters serve to demarcate a region of uninterrupted growth from that of collapse [16].

In the current study, we consider a model assembly system that shows the unstable dynamics of assembly around the critical concentration of ATP/GTP. The present model utilizes CMT for describing the behavior of monomers in the solvent and polymer for simplification of analysis. We applied a kinetic model that unifies de novo nucleation and growth by considering the monomer-monomer interactions as a diffusion process. In addition, the diffusion process of the monomeric protein has been considered from the perspective of non-linearity. According to Fick’s law, the continuity of monomer concentration of c (i = 1, n) including chemical reaction items, may be described using diffusion coefficients D , kinetic coefficients k , and the concentrations of individual compounds c . Protein assembly is limited by the slow diffusion rate of monomer proteins, which is a diffusion-rate limiting aggregation process. Therefore, diffusion items and reaction items cannot necessarily be separated; therefore, we described kinetic rate of c as follows [8]:

Here, the first item on the right represents the diffusion rate. The second item, f (c ), denotes the function of kinetic rate of reactions other than the diffusion process. k represents a coefficient.

Methods

Numerical simulation Numerical calculations were performed using Mathematica 8 (Wolfram Research, Inc., Champaign, IL).

Results

General formulation of an assembly

The model consists of several steps: (i) the monomer achieves an interactive state by binding a cofactor (ATP/GTP) that provides the monomer with the ability to interact; (ii) the monomer itself possesses the ability to hydrolyze the cofactor and lose assembly activity; (iii) the monomer has the ability to exchange the inactive hydrolyzed cofactor (ADP/GDP) with an active non-hydrolyzed one; and (iv) ATP/GTP are supplied continuously from the external environment. The second requirement indicates a self-limiting property of the monomer that causes dynamic instability during monomer-monomer interaction. When examining protein interaction kinetics, analysis of the fluctuation in monomer concentrations was performed using Mathematica 9.

Protein interaction kinetics

The model scheme is shown in Fig. 1. There are three types of monomer: ATP/GTP-binding monomer X, ADP/GDP-binding monomer Y in the oligomer (W), and the released ADP/GDP-binding monomer Z. X has the higher assembly activity, and Y and Z have lower assembly activity. We set the oligomer concentration W to be a constant, as de novo assembly is considered much slower than monomer interaction in the steady state [11-14]. The individual steps are shown below:

Fig. 1

Scheme of monomer interaction. Individual globules or oblongs represent monomers X, Y, Z, and oligomer W. Kinetic coefficients, k , k , k , and k are shown next to the arrows. Outside and inside signify the outside and inside of the cell, respectively. Y is located at the end of the oligomer W

First, X associates with the assembly nucleus W to be Y at the end of W.

In the next step, the intermediate species Y is released to be Z:

Z recovers its interaction activity by exchanging the active cofactor ATP/GTP (P) for the inactive cofactor ADP/GDP (P’), returning to X (see Fig. 1):

In addition, direct slow conversion is supposed:

The kinetic equations were set according to the simple reaction cascade described above. We obtained equations for the protein interaction kinetics using the diffusion coefficient:

In addition, the total concentration of the monomer is maintained constant.

M, which represents the total concentration of the monomeric proteins, is maintained constant. A simple consideration of the diffusion-limited step implies that, when the kinetic rate can be described according to Fick’s law using the diffusion coefficients D , D and D then [17-19]:

As the oligomer diffusion rate is small, we set D , D > >D . Therefore, m , and m are substantially proportional to D and D , respectively. Accordingly. kinetic coefficients k and k were defined as the proportional coefficients below:

Rewriting (6), (7), and (8) using (10) and (11),

In the above equations, k and k represent the kinetic coefficients for the addition of the monomer to the oligomer and the release of the monomer from the oligomer, respectively.

In order to obtain the monomer concentration at the steady state of the reaction system, the right-hand side of Eqs. (12), (13), and (14) were set to be equal to zero and Eq. (9) were used to give:

In the above approximation, we omitted D D and k as the diffusion coefficients and the direct conversion rate of X into Z is small.

Fluctuation of diffusion coefficient

Next, we considered the fluctuations of participant proteins using small letters x, y, and z:

In Eq. (14), the subscript ‘e’ signifies values at the steady state.

In an assembly, monomers associate with other monomers. From Eq. (9),

Therefore, the fluctuation y may be represented using −x−z. The fluctuation kinetics are thus provided by two parameters, namely x and z.

Given the nonlinearity during diffusion, we assume kinetic instability in the monomer-monomer interaction, and that the sensitivity of the assembly in response to environmental change may be evaluated. Indeed, the diffusion coefficient D of one macromolecule in the solution may generally be represented using the fluctuation concentration c :

c denotes the concentration of the solute, α is a coefficient, and D is the diffusion coefficient when the fluctuation of monomeric protein is negligible. The dependence of the diffusion coefficient on the protein concentration has been reported [20, 21]. O’Learly reported that diffusion coefficients of proteins linearly decrease in proportion to the concentration, when the latter is sufficiently small. The fluctuation of the diffusion coefficient is obtained by considering the dependence of the coefficients on the concentration of the monomer from Eq. (20) [8]:

Here, the fluctuation term αx (α >0) and γx (γ >0) contributes to a decrease in D and D , as higher assembly activity reduces diffusion. In contrast, an increase in the fluctuation terms βz (β >0) and δz (δ >0 serves to increase the diffusion coefficients D and D , as lower interaction or assembly activity increases diffusion. When the assembly activity of Z is lower, the fluctuation item δz is negligible, in accordance with the fluctuation kinetic equations given by (19), Eqs. (12), (14), (21), and (22):

Here, y, fluctuation of intermediate species Y is negligible as the value is sufficiently small. In addition, we used the following equations to describe the balance in detail:

and

To simplify the notation in (23) and (24), we set:

and obtained:

Eqs. (28) and (29) represent a master equation for the application of CMT.

Calculus simulation of concentration oscillations

For analysis of the behavior of the system, including multi-parameters, the examination of the linearization of behavior of the system near a steady state provides insights into the qualitative behavior of the system in the neighborhood of the point. In particular, the eigenvalues of the linear part of the governing kinetic equations enable determination of the stability of the system behavior. CMT is a rigorous formulation of this observation that enables the reduction of a large number of parameters [22].

Around the steady state (x, z) = (0, 0) of Eqs. (28) and (29), the Jacobian matrix of (dx/dt, dz/dt) is given by:

Subsequently, the time-course of the monomer concentrations was simulated by substituting appropriate numerical values into Eqs. (28) and (29). The simulation results under the above conditions are shown in Fig. 2. A numerical calculation was performed over a sufficiently long period to evaluate the assembly trend. The steady-state concentrations of X and Z are given by Eqs. (15) and (17). The critical value of p is given by

Fig. 2

Time-course of the fluctuation in monomer concentrations displays a oscillation. Diffusion of active cofactor binding monomer (X) and of inactive cofactor binding monomer (Z). p is (a) 0.000, (b) 0.001, (c) 0.002, (d) 0.004, (e) 0.008, (f) 0.009, (g) 0.01000, (h) 0.010705, (i) 0.011000. The graphs show plots of X (black), Y(red), and Z (blue). Lines represent the concentration of X and Z. The horizontal axis represents time (0 ≤ t ≤ 1000) and the vertical axis represents the concentration of X and Z. When p exceeds 0.01, oscillations are observed. The Mathematica (version 9, Wolfram Research, Inc., Champaign, IL) code for p = 0.01 is as follows: p = 0.01 X = ((D2 M p)/(D2 k + D2 p + D1 D2 W + D1 p W)) Y = ((D1 M p W)/(D2 k + D2 p + D1 D2 W + D1 p W)) Z = ((D2 M (k + D1 W))/(D2 k + D2 p + D1 D2 W + D1 p W)) M = 0.1 W = 1 D1 = 0.28 D2 = 0.012061855670103093` a = 150 b = 156 k = 0.005 c = 0.1 d = 0 NDSolve[{Derivative[1][x][t] == − (D1 - a X) x[t] + a x[t]^2 + (p - b X) z[t] - b x[t] z[t] - k x[t], Derivative[1][z][t] == k x[t] + c x[t]^2 + d x[t] z[t] - p z[t], x[0] == 1.`*^-6, z[0] == 1.`*^-6}, {x, z}, {t, 0, 3300}, MaxSteps - >50000] g001 = Plot[{X + x[t]} /. %, {t, 0, 1000}, PlotRange - > All, PlotStyle - > {RGBColor[0, 0, 0]}] g002 = Plot[{Y - x[t] - z[t]} /. %%, {t, 0, 1000}, PlotRange - > All, PlotStyle - > {RGBColor[1, 0, 0]}] g003 = Plot[{Z + z[t]} /. %%%, {t, 0, 1000}, PlotRange - > All, PlotStyle - > {RGBColor[0, 0, 1]}, PlotRange - > All] Show[g001, g002, g003]

Here, the small affix c indicates the critical point of ATP/GTP concentration. Next, we conducted a simulation with values of M = 0.1, X  = 0.002, D  = 0.28, D  = 0.012, a = 150, b = 150, k = 0.005, c = 0.0, and d = 0. Solving the above with respect to p with substitution of these values in Eq. (31), we find the critical value:

As a result, the fluctuations oscillate between decrease and increase in monomer concentrations, as shown in Fig. 2. When p