Kinetic Analysis of Protein Assembly on the Basis of the Center Manifold around the Critical Point
1 Kinetic Analysis of Protein Assembly on the Basis of the Center Manifold around the Critical Point
Tatsuaki Tsuruyama* *
Department of Diagnostic Pathology, Kyoto University Graduate School of Medicine 1 Yoshida-Konoe-cho, Sakyo-ku, Kyoto, Kyoto Prefecture 606-8501, Japan Phone: +81-75-753-4427; Fax: +81-75-753-4493 E-mail: [email protected] Abstract
Protein assembly plays an important role in the regulation of biological systems. The cytoskeleton assembly activity is provided by the binding cofactors GTP (guanidine triphosphate) or ATP(adenosine triphosphate) to monomeric protein, and is initiated by assembling the monomeric proteins. The binding GTP or ATP is hydrolyzed to GDP (guanidine diphosphate) or ADP (adenosine diphosphate) by the monomeric enzymatic activity. This self-limited assembly is characteristic of the cytoskeleton. To quantitatively evaluate the assembly kinetics, we propose a nonlinear and non-equilibrium kinetic model, with the nonlinearity provided by the fluctuation in monomer concentrations during the diffusion. Numerical simulations suggest that the assembly and disassembly oscillates in a chaos-like manner. We use a kinetic analysis of the center manifold around the critical point to show that minimal increases in ATP/GTP concentrations may lead to some attenuation in the amplification of these fluctuations. The present model and our application of center manifold theory illustrate a unique feature of protein assemblies, and our stability analysis provides an analytical methodology for the biological reaction system.
Keywords: protein assembly, nonlinear kinetics, fluctuation
1. Introduction
Protein assembly is essential for cellular activities such as cytoskeleton formation and protein assembly [1, 2] . A transient binding cofactor such as ATP/GTP to monomeric protein controls the assembly of proteins such as actin and tubulin. Various mathematical models have been proposed to explain the kinetic signaling behavior based on a set of kinetic equations [3]. In the current study, we aim to understand the process of protein assembly using the following two novel perspectives, which are formulated on the basis of a nonlinear and non-equilibrium kinetic model. In fact, cytoskeletal protein assembly is a dynamic process between assembly and disassembly [4-8]. According to Fick’s law, the continuity of the monomer concentration of c i ( i = 1 ,…, n ) including chemical reaction items can be described using diffusion coefficients D i , kinetic coefficients k i , and the concentrations of individual compounds c i . Protein assembly is limited by the slow diffusion rate of monomer proteins, which is a diffusion-limited aggregation process. Therefore, diffusion items and reaction items cannot necessarily be separated and we therefore described kinetic coefficients that depend on the diffusion coefficients as follows[8]: dc i dt = k i D i c i + k i c ij = n ∑ c j + ! (1) We consider a model system that is far from equilibrium because of a continuous supply of chemicals from the outside. The chemical species are GTP and ATP, cofactors that provide monomers with assembly activity. The aim of this study is to illustrate the behavior of non-equilibrium protein assembly. Materials and methods
The model consists of several steps: (i) the monomer achieves an interactive state by reversibly binding a cofactor (ATP/GTP) that provides the monomer with interaction activity; (ii) the monomer itself has the ability to hydrolyze the cofactor and lose the assembly activity; (iii) the monomer has the ability to exchange the inactive hydrolyzed cofactor (ADP/GDP) with an active non-hydrolyzed one; (v) ATP/GTP are supplied continuously from the outside. The second requirement indicates one of the self-limiting properties of the monomer, causing dynamic instability in the monomer interaction. In the protein interaction kinetics, our analysis of the fluctuation in monomer concentrations was performed using Mathematica 9 (Wolfram Research, Inc., Champaign, IL). The model scheme is shown in Figure 1. There are three types of monomer, ATP/GTP -binding monomer X , ADP/GDP--binding monomer Y in the oligomer ( W ), and released ADP/GDP-binding monomer Z . X has the higher assembly activity, and Y and Z have the lower assembly activity. We set the oligomer concentration W to be a constant, because de novo assembly is considered to be much slower than monomer interaction in the steady state [6, 7, 12-15]. First, X associates with the assembly nucleus W to be Y at the end of W . X + W → W+Y ( m ; kinetic coefficients) (2) In the next step, the intermediate species Y is release to be Z : Y → Z ( m ) (3) Z recovers its interaction activity by exchanging the active cofactor ATP/GTP ( P ) for the inactive cofactor ADP/GDP ( P ’), returning to X (see Figure 1): Z+P → X+P’ ( k ) (4) In addition, direct slow conversion is supposed: X → Z ( k ) (5) The kinetic equations were set according to the simple reaction cascade described above. We obtain equations for the protein interaction kinetics using the diffusion coefficient: dX/dt = − m WX + k PZ − k X (6’) dY/dt = m WX – m Y (7’) dZ/dt = m Y − k PZ + k X (8’) X + Y + Z = M (9) 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 X and D Y , then[9-11]: m ∝ D X + D W = k D (10) m ∝ D W + D Y = k D (11) Rewriting (6’)-(8’) using (10) and (11), dX/dt = − k D WX + k PZ − k X (6) dY/dt = k D WX – k D Y (7) dZ/dt = k D Y − k PZ + k X (8) To obtain the monomer concentration at the steady state of the reaction system, setting right-hand side of Eqs. (6) – (8) equal to zero and using Eqs. (9) gives: X e = D k MpD k k + D k p + D D k k W + D k pW ~ D k MD k + D k W (12) Y e = D k MpWD k k + D k p + D D k k W + D k pW ~ D k MWD k + D k W (13) Z e = D k M k + D k W ( ) D k k + D k p + D D k k W + D k pW ~ D k M k + D k W ( ) D k p + D k pW ~ 0 (14) Here, the item D D and k are negligible because the diffusion coefficients are small and the direct conversion rate of X into Z is small. We now consider the fluctuation of participant proteins: X = X e + x , Y = Y e + y, Z = Z e + z (15) In Eq. (14), the subscript ‘e’ signifies values at the steady state. In an actual assembly, monomers associate with other monomers. From Eq. (9), we have x + y + z = 0 and therefore the fluctuation y can be written using – x–z. Therefore, the fluctuation kinetics are provided by two parameters, X and Z . Given the nonlinearity during diffusion, we assume there is kinetic instability in the monomer interaction, and that the sensitivity of the assembly in response to the environmental change can be evaluated. Indeed, the diffusion coefficient D of one macromolecule in the solution can generally be written as: D = kT/ η (1 – α c) ≈ kT (1 + α c) / η = D (1 + α c) (16) where T is the temperature of the solution, k is the Boltzmann constant, η is the frictional coefficient of a given macromolecule in solution and c denotes the concentration of the solute, and α is a coefficient, and D is the diffusion coefficient [16-20] . The fluctuation of the diffusion coefficient is given by considering the dependency of the coefficients upon the concentration of the monomer [8]: dD x = α x – β z (17) dD z = γ x – δ z (18) Here, an increase in X contributes to a decrease in D x and D z in the fluctuation term α x ( α > 0) and γ x ( γ > 0), because the higher assembly activity reduces diffusion. In contrast, an increase in Z serves to increase the diffusion coefficients D x in the fluctuation terms β z ( β > 0) and δ z ( δ > 0), because the lower interaction or assembly activity increases diffusion. In fact the assembly activity of Z is lower, the fluctuation item δ z is negligible. Because x + y + z = 0, Eqs (6), (8), (17), and (18) give the fluctuation kinetic equations : dx/dt = ( − k D x W+k α WX e − k ) x +( − k β WX e + k P ) z + k W α x − k W β xz (19) dz/dt = ( k − k γ Y e ) x − k Pz + k γ x + k γ xz (20) To simplify the notation, we set: k D x = D , k β = b, k α = a, k γ = c , k P = p, k = k (21) and obtain: dx/dt = ( − D W+aWX e − k ) x + ( − bWX e + p ) z + Wax − Wbxz (22) dz/dt = ( k − cY e ) x − pz + cx + cxz (23) Around the steady state ( x , z ) = (0, 0) of Eqs. (22) and (23), the Jacobian matrix of ( dx/dt , dy/dt ) is given by: L = − D W + aX e − k − bWX e + pk − cY e − p ⎡⎣⎢⎢ ⎤⎦⎥⎥ (24)
3. Results 3.1 Calculus simulation of concentration oscillations
The time-course of the monomer concentrations was simulated by substituting appropriate numerical values into Eqs. (22) and (23). 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. (12) and (14). The simulation results under the above conditions are shown in Figure 2. When p varies about the value given by det [ L ]= ( − D W+aX e − k )( − p ) − ( k − cY e ) ( − bWX e + p ) = 0 (25) Next, we conducted a simulation with values of M = 0.1, W = 1, D = 0.28, D = 0.012, a = 150, b = 150, k = 0.005, c = 0.1, and d = 0. Solving (24) with respect to p with substitution of these values into (24), we find p ≈ p c (26) The small affix c indicates the critical point of ATP/GTP concentration. As a result, the fluctuations oscillate between decreasing and increasing as shown in Figure 2. When p < p c , the fluctuation was found to attenuate (Figure 2 d ) and the monomer concentration reached the plateau. While, when p > p c , the fluctuation was found to diverge (Figure 2 f ). We considered the kinetic stability around the critical point ( p = p c ) using minimal difference ε . The matrix L is then given by: Lc = − D W + aX e − k − bWX c + p c k − cY e − p c ⎡⎣⎢⎢ ⎤⎦⎥⎥ (27) Using the eigenvectors of Lc , [ l l ], we perform the following coordinate transformation: ! x ! z ⎡⎣⎢ ⎤⎦⎥ = l l ⎡⎣ ⎤⎦ ! u ! v ⎡⎣⎢ ⎤⎦⎥ (28) ! u ! v ⎡⎣⎢ ⎤⎦⎥ = l l ⎡⎣ ⎤⎦ − ! x ! z ⎡⎣⎢ ⎤⎦⎥ (29) In reference to the numerical simulation (Figure 2), when D , k, p are sufficiently small, [ l l ] ≈ [{-( aWX/k ), 1}, { aWX/k , 1}]. We subsequently set: du / dt = f u ( u, v ) (30) dv / dt = f v ( u, v ) (31) The center manifold around the critical point ( p = p c ) is then given as follows : u= h ( ε , v ) = a v + a v ε + a ε + a v + a v ε + a v ε + a ε + O ( ε ) (32) We then analyzed the effect of changes ε in p using the center manifold around the critical point of the system. Subsequently, u = ( dv/dt ) ∂ h ( u, ε )/ ∂ u + ( d ε / dt ) ∂ h ( u, ε )/ ∂ ε = (2 a v + a ε ) f u ( u, v ) (33) Using Eqs. (31) and (32), we then obtain: (2 a v + a ε ) f u ( u,v ) = a v + a v ε + a ε + a v + a v ε + a v ε + a ε + O ( ε ) (34) Solving Eq. (34) gives the coefficients of a i in Eq. (32): a = a = 0. Substituting u in Eq. (32) given by ν and ε into f v ( u,v ) in Eq. (29), we can obtain the kinetic stability equation for fluctuation ν using the coefficients n i ( i = 1, …, 7) as follows: dv/dt = n v + n v ε + n ε + n v + n v ε + n v ε + n ε + O ( ε ) (35) Independent of the numerical values in Eq. (35), n , n , n = 0 (36) Then, we have: dv/dt = n v + n v ε + n v + n v ε + O ( ε ) (37) By setting left-hand side equivalent to zero, v = − n − n ε ± ( n + n ε ) − n n ε n (38) We have an approximate solution to Eq. (38): v = 0, ≈ c ε , − n /n (39) Here, c is a constant coefficient and two of the oscillatory amplitude of v except – n / n are nearly equivalent to zero when the concentration of ATP/GTP is equivalent to critical concentration ( ε ≈ c ε increases as the concentration of ATP/GTP increases. From (32), we have formulation of u using a constant coefficient c’ , u ≈ c’ ( n /n ) (40) When D , k, p are sufficiently small, substituting [ l l ] = [{ − aWX/k , 1}, { aWX/k , 1}] into (27) gives : x = − ( aWX /k ) u + v ≈ v (41) As a result, as we described v, two of the amplitudes of x are nearly equivalent to zero ( ε ≈ − n / n ) increases as ε increases (Figure 2). Thus, stability analysis can predict the behavior of the fluctuation around the critical point of protein assembling system.
4. Discussion
We have presented a model for the protein assembly kinetics, and analyzed the stability around the critical point using center manifold theory. The nonlinear kinetic equations include three parameters (
X, Y, and Z ), but only two are actually independent (Eq.(8)). In the simulations, ATP/GTP or ADP/GDP binding monomers periodically exhibit a chaotic oscillation between assembly and disassembly. This accurately reflects the microtubule kinetics showing instable assembly [4, 5]. To the best of our knowledge, this is one of the first reports on the application of center manifold theory to the analysis of biological reaction systems [8]. The fluctuation of monomer concentrations was subjected to a perturbation expansion using a minimal increase in the supply of ATP/GTP near the concentration at the critical point. This mathematical method can precisely treat nonlinear and multi-parameter systems around the critical point. The fluctuation kinetics are expected to change from convergence to divergence of the concentration fluctuation of the monomer, i.e., stable to instable, around the critical point as shown in Figure 2. Because of this high sensitivity to the concentration of ATP/GTP, the protein assembly is dynamically regulated by minimal changes in the supply of ATP/GTP, which is controlled metabolically. Still, the theoretical basis of protein assembling for quantitative evaluation requires further formulation on the basis of experimental study. Acknowledgement
This work was supported by a Grant-in-Aid for Scientific Research on Innovative Areas, MEXT, Japan. Project name is ‘‘Synergy of fluctuation and Structure’’. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. References
1. Arnal I, Karsenti E, Hyman AA (2000) Structural transitions at microtubule ends correlate with their dynamic properties in Xenopus egg extracts. J Cell Biol 149: 767-774. 2. Woolley DM (2010) Flagellar oscillation: a commentary on proposed mechanisms. Biol Rev Camb Philos Soc 85: 453-470. 3. Hazra P, Inoue K, Laan W, Hellingwerf KJ, Terazima M (2006) Tetramer formation kinetics in the signaling state of AppA monitored by time-resolved diffusion. Biophys J 91: 654-661. 4. Wu Z, Wang HW, Mu W, Ouyang Z, Nogales E, et al. (2009) Simulations of tubulin sheet polymers as possible structural intermediates in microtubule assembly. PLoS One 4: e7291. 5. VanBuren V, Cassimeris L, Odde DJ (2005) Mechanochemical model of microtubule structure and self-assembly kinetics. Biophys J 89: 2911-2926. 6. Symmons MF, Martin SR, Bayley PM (1996) Dynamic properties of nucleated microtubules: GTP utilisation in the subcritical concentration regime. J Cell Sci 109(11): 2755-2766. 7. Voter WA, Erickson HP (1984) The kinetics of microtubule assembly. Evidence for a two-stage nucleation mechanism. J Biol Chem 259: 10430-10438. 8. Tsuruyama T (2014) A model of cell biological signaling predicts a phase transition of signaling and provides mathematical formulae. PLoS One 9: e102911. 9. Wustner D, Solanko LM, Lund FW, Sage D, Schroll HJ, Lomholt, MA (2012) Quantitative fluorescence loss in photobleaching for analysis of protein transport and aggregation. BMC Bioinformatics 13: 296. 10. Dorsaz N, De Michele C, Piazza F, Foffi G (2010) Inertial effects in diffusion-limited reactions. J Phys Condens Matter 22: 104116. 11. Kasche V, de Boer M, Lazo C, Gad M (2003) Direct observation of intraparticle equilibration and the rate-limiting step in adsorption of proteins in chromatographic adsorbents with confocal laser scanning microscopy. J Chromatogr B Analyt Technol Biomed Life Sci 790: 115-129. 12. Carlier MF (1990) Actin polymerization and ATP hydrolysis. Adv Biophys 26: 51-73. 13. Carlsson AE (2008) Model of reduction of actin polymerization forces by ATP hydrolysis. Phys Biol 5: 036002. 14. Korn ED, Carlier MF, Pantaloni D (1987) Actin polymerization and ATP hydrolysis. Science 238: 638-644. 15. Pantaloni D, Hill TL, Carlier MF, Korn ED (1985) A model for actin polymerization and the kinetic effects of ATP hydrolysis. Proc Natl Acad Sci USA 82: 7207-7211.
16. Carles H, Ziya A (1981) Concentration dependence of diffusion coefficient at various molecular weights and temperatures. Polymer 22: 1165–1168. 17. Yamahawa H (1971) Modern Theory of Polymer Solutions. New York. Harper& Row, 266–267. 18. Gallagher WH, Woodward CK (1989) The concentration dependence of the diffusion coefficient for bovine pancreatic trypsin inhibitor: a dynamic light scattering study of a small protein. Biopolymers 28: 2001–2024.
19. Nesmelova IV, Skirda VD, Fedotov VD (2002) Generalized concentration dependence of globular protein self-diffusion coefficients in aqueous solutions. Biopolymers 63: 132–140. 20. O’Leary TJ (1987) Concentration dependence of protein diffusion. Biophys J 52: 137–139. Figure Legends
Figure 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. Figure 2. Time-course of the fluctuation in monomer concentrations displays a chaotic 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.010, (g) 0.010705, and (h) 0.011. The upper graph shows two parametric plots of X (red), and Z (blue). Lines represent the concentration of X and Z . The horizontal axis represents time (0 ≤ t ≤ X and Z . When p exceeds 0.01, chaotic 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, Figure 1 Figure 219