Open access peer-reviewed chapter

Non-Linear Kinetic Analysis of Protein Assembly Based on Center Manifold Theory

By Tatsuaki Tsuruyama

Submitted: April 18th 2017Reviewed: August 30th 2017Published: December 20th 2017

DOI: 10.5772/intechopen.70750

Downloaded: 708


This review introduces a novel mathematical description of protein assembly. Protein assembly occurs in a substantially open non-equilibrium and non-linear kinetic system. The goal of systems biology is to make predictions about such complicated systems, but few have conducted stability analysis for given systems. Particularly, simulated dynamic behaviors have not been sufficiently verified by kinetic analysis in predicting macromolecular protein interactions and assembly. The non-linearity of protein assembly kinetics is complex, and it is very difficult to determine a model of multi-protein interactions based on numerical calculation. We studied the non-linear kinetics involved in the diffusion process of proteins consisting of two or three species of macromolecules and set a novel model in which non-linearity is given by the diffusion coefficient that depends on the protein concentration. By making the diffusion coefficient concentration-dependent, non-linearity leads to a simple system model. Protein assembly is initiated by monomeric protein interactions and regulated by cofactors such as guanidine triphosphate (GTP) or adenosine triphosphate (ATP) binding to the monomer. This cofactor concentration promotes the dynamic behavior of protein assembly and can be treated as an order parameter. Further, kinetic stability analysis in the center manifold theory (CMT) is introduced for analyzing the behavior of the system around the critical state. Although CMT has not been sufficiently applied for stability analysis of protein assembly systems, this theory predicts the dynamic behavior of the assembly system around the critical point using concentration as a cofactor. Protein assembly theory will provide a novel framework for nonlinear multi-parametric analysis.


  • protein assembly
  • center manifold theory
  • tubulin
  • non-linear kinetics
  • non-equilibrium state
  • diffusion coefficient
  • oscillation

1. Background

Protein assembly is essential for cellular activities such as cell signaling, gene expression by transcription factor complexes, cytoskeleton formation, endocytosis, and cell motility. This reaction system is one of the complicated systems and non-linear kinetics has been applied for understanding the dynamic behavior [1, 2, 3, 4, 5, 6, 7, 8]. Among the protein assembly, tubulin and actin polymerization are well-known events that have been analyzed using numerous methods [9, 10, 11, 12, 13]. Microtubule and actin filaments consist of monomers that bind the cofactors guanidine triphosphate (GTP) or adenosine triphosphate (ATP) to acquire assembly activity [14, 15, 16]. In general, a protein interaction is controlled by an electric charge on the amino acid residue(s) of the protein, such as tyrosine, serine, and threonine, by covalently binding to or reacting with a cofactor such as ATP or GTP; subsequently, the monomer loses its interaction activity by hydrolyzing ATP or GTP into ADP or GDP or through dephosphorylation, which is mediated by the other protein’s phosphatase activity or its own enzymatic activity [17, 18].

In particular, dynamic instabilityin tubulin polymerization has been extensively investigated. Dynamic instabilitysignifies the intermittent transition between slow growth and rapid shrinkage in polymeric assemblies of microtubules [9, 10, 11, 12, 13, 19]. Further, intra-polymeric Brownian motion and fluctuation influence the structure and elasticity of tubules [20]. Zapperi and Mahadevan presented an excellent model in which the ratio of longitudinal to lateral interactions characterizes the assembly [21]. Hammele et al. presented a physical model and suggested the physical properties of the microtubules [22]. Nucleation is the rate-limiting step controlling the overall polymerization process [18]. The stable nucleus for polymerization is oligomers, and the growth of aggregates through elongation/dissociation follows. For stable growth, tubule lifespan is controlled by a GTP-cap that forms at their ends [19].

As another example of protein interactions, in the mitogen-activated protein kinase (MAPK) signaling cascade, a set of protein kinases and protein substrates construct the signaling network. The cofactor ATP/GTP transfers biological information in the reaction network to alter gene expression [6, 7, 8]. Mathematical models of this cascade have demonstrated that the system can act as an ultra-sensitive switch based on a combination of phosphorylation of protein substrates and implicit feedback, leading to multi-stability [23, 24, 25]. Recently, Ueno et al. reported that a model of MAPK signaling cascade functions as a band-path filtering system.


2. Protein interaction kinetics

2.1. Protein interaction model

Steps in protein assembly of microtubular polymerization are summarized as follows: (i) the protein achieves an interaction active state by reversibly binding to a cofactor, which provides the protein with assembly or interaction activity; the protein interaction activity decreases when a hydrolyzed inactive cofactor is bound compared to an active cofactor; (ii) the protein can hydrolyze the cofactor; (iii) the protein can exchange the inactive cofactor, such as ADP, with an active cofactor; and (iv) active cofactors are supplied continuously and externally. Thus, the interaction activity is self-limiting, in which the protein itself limits assembly activity, resulting in dynamic instability [10, 11, 12, 13, 19, 23, 26, 27, 28].

Let us consider a three-monomer model in which an active cofactor-binding protein (X), oligomeric protein (W), and inactive cofactor-binding protein (Z) coexist. The sum of monomeric proteins is kept constant; this kinetics is attributed to a two-parametric analysis. First, Xcan reversibly associate with the oligomer:

X+WWkinetic rate coefficientsk1E1

Subsequently, Zis released from the oligomer:


Here, Pi signifies the released inorganic phosphate. Finally, in the presence of sufficient active cofactor, Zreleases the inactive cofactor, ADP or GDP (I), binds an active cofactor, ATP or GTP (P), and recovers its interaction activity, returning to the protein


In addition, a monomeric protein has the potential to hydrolyze the cofactor by interaction:


Formula (5) represents a self-reproducing reaction, which yields non-linear kinetics in this reaction system.

2.2. Concentration dependence of protein diffusion

Diffusion of proteins plays an important role in protein interactions. Analysis of dilute solutions of a macromolecule requires a greater understanding of the concentration dependence of the diffusion rate because of the hydrodynamics of protein solutions involving mutual diffusion of protein molecules. One of the approaches is applicable to the study of self-diffusion in solutions [29, 30, 31, 32]. In fact, proteins interact or associate with other monomeric proteins and phosphorylate or are phosphorylated by the proteins. In dilute solution, proteins may diffuse in a free manner with sufficiently large vacant space that accounts for only a fraction of the volume of a protein molecule. These vacancies are sufficiently large to be occupied by proteins that are as large as the hydrodynamic volume. The effects of molecular shape and size, solvent, and environment such as ion intensity and pH determine the concentration dependence of the diffusion rate. Here, the diffusion rate D0 is set as the probability of vacancy formation and does not depend on the velocity at which a monomer diffuses. The probability Pis given by a void adjacent to the objective protein, which is sufficiently large to permit diffusion: D = DoP.

The probability P(V) of forming a vacancy of volume V in the solution is given by


where Veis protein exclusion volume, c is the protein solution concentration, and β is a constant that reflects the effects of interactions with other macromolecules and shape on the probability of vacancy formation. Thus, Ddepends on the size and shape of the protein and likely other factors


At low concentrations of the macromolecule


Thus, the dependence of the diffusion coefficient can be described using the protein solution concentration (Figure 1).

Figure 1.

Self-diffusion rate constantsDfor concentration of protein,c, with theoretical curve fromEq. (8)with β = 3.0. Near the point x ~ 0, the diffusion coefficient obeysD = 1 − 2c.The vertical axis represents the diffusion coefficientD × 107 (cm2/s) and horizontal axis represents hemoglobin concentration (g/dL). The graph is shown using new arbitrary values with reference to experimental data.

2.3. Viscosity and diffusion coefficient of protein

The compatibility of Eq. (8) with the experimental data strongly suggests that the concentration dependence of the protein diffusion constant is governed by excluded volume interactions, which may be predicted by calculating equilibrium protein density fluctuations. Eq. (8) is consistent with the equation describing the viscosity ηof concentrated protein solutions


where η0 is the solution viscosity at infinite dilution, η is the intrinsic viscosity of the solution, cis the protein solution density, and k/vis a constant that corrects for the overlap of free volume v.

2.4. Diffusion-limited protein association

Because protein diffusion in the viscous cytoplasm is significantly slower than other reactions between interactions with other molecules, the protein assembly reaction is a diffusion-controlled or diffusion-limit process. In general, the diffusion rate of monomer is given by the concentration gradient that is the molar flux per unit area, J. This is the additional rate of monomer Xtowards oligomer Wmultiplied by the area of the spherical surface of radius Rof the reactive end of polymer or oligomer

Rate of reaction=4πR2JE10

In Fick’s first law, the flux towards Xis proportional to the concentration gradient at radius of macromolecule R


By substitution of Eq. (10) into Eq. (11), the rate of reaction v is given as


The rate of the diffusion-controlled reaction is equal to the average flow of Xmolecule to all Wmolecules. Accordingly, the global flow of all Xto Wis 4π R*DXNA XW.Similarly, the flow of all Wto Xis 4π RDWNA WX.Further, using the sum of the diffusion coefficients of the two species, the diffusion coefficient is rewritten as D = (DX + DW)/2. Then, the addition rate of Xto Wis given by

Addition rate=4πRkXoDXWE13

In reality, the diffusion rate of oligomer is negligible relative to that of monomer, and the addition rate is given as

Addition rate=4πRkXoDXXW=kXDXXWE14



Using the above formula, the kinetic rate of ciand items of interaction between monomers is as follows [8]:


where f(cj) represents the sum of the reaction kinetic items of jth species, cj, except the diffusion process. When the potential of the electric field, ϕ(r), is considered, the flux is described using spatial coordinate r


In the solvent including two chemical species, the relative movement is related to the local concentration gradient of cjand potential energy. When the gradient is described using Eq. (17), monomer Xmoves across the sphere with radius Rsurrounding Wand can be described by


when Xand the active site on the oligomer or polymer Winteract to assemble or elongate the polymer, which is determined by R.The total flux will be equivalent to the chemical reaction rate using an arbitrary kinetic coefficient k


At the steady state, flux across the sphere with a radius, or a shape parameter, ris constant for any values. Accordingly, Rin Eq. (19) is replaced with rusing Eq. (18)


By integration and rearrangement of Eq. (20)


and here


Because r- > ∞, V(r) approaches zero


Here, Xr=R = (k/kR) X.In calculating the above formula, fluctuations of diffusion coefficients are neglected when the concentration of Xis sufficiently low and kept constant during the reaction. The kinetic coefficient is controlled by the diffusion process and the aggregation or assembly is completed through interactions between Xand W; Xr=R is set to zero. Then, kis given by Eq. (23):


The diffusion coefficients can be altered in proportion to the fluctuation of monomeric protein concentration [9, 10, 11, 12, 13, 14]. In the derivation of Eq. (23), the diffusion coefficient is related to the viscosity ηby the Einstein-Stokes formula


By using the Gibbs-Duhem expression, the diffusion coefficient Dof one macromolecule can be written as


where Tis the temperature of the solution, kBis the Boltzmann constant, and η1 is the frictional coefficient of the macromolecule in solution. A1 is the second virial coefficient, v1 is the partial specific volume of protein with molecular weight M1, and NAis Avogadro’s number. The small letter c1 denotes the concentration of the solute macromolecule. Then, dependency of the diffusion coefficient on the ith component, Di, is as follows from (26):


where vjis the partial specific volume of the polymer with molecular weight Mj.

Further, the diffusion coefficient is given by extending the above formula to a mixed solution of two macromolecules, Xand Z


where vXand vZare the partial specific volumes of Xand Zwith molecular weights MXand MZ, respectively. AXand AZare the second virial coefficients.

2.5. Fluctuation of diffusion coefficient

Subsequently, let us consider the fluctuation of participant proteins using lowercase letters


Here, the sum of the kinetic rate of fluctuations is constant because total amount of monomeric protein is constant


The dependency of the diffusion rate refers to the high interaction activities of Xand low interaction activities of Z. An increase in Xcontributes to a decrease in the diffusion coefficients DXand DZin fluctuation items because of the higher interaction activity that reduces diffusion; in contrast, an increased Zcontributes to increased diffusion coefficients because the interaction between increased Zinduces lower assembly activity for the interaction with monomeric proteins, resulting in increased mobility of monomeric proteins. This dependency gives the non-linearity fluctuation items in the kinetic equation. Here, all preparations for kinetics are completed.

2.6. General theory of non-linear kinetic equation of protein assembly

According to the above simple reaction cascade, (1)-(5), the kinetic equation contains the concentrations of monomeric proteins as variables. For simplification, the equations are written as follows:




Further, for simplicity, Eqs. (34) and (36) are given by replacing the kinetic coefficients with arbitrary coefficients


Here, k1D1W = k1, k4D4W = k4, k5D5 = k5, and p = k3 P.At the steady state, setting the right hands of Eqs. (32)-(34) equal to zero


Here, the fluctuation of diffusion coefficients are given by


By altering X, Z, k1, k4’,and k5into X + x, Z + z, k1’ – ax + bz, k4’ − cx + dz, and k5’ − ex + fzin Eqs. (36) and (37) and arranging, two dependent equations are obtained


Accordingly, the overall behavior of the kinetics of protein assembly is given by the monomeric kinetics of xand z. fxz(), fzx(), fxx(), and fzz() represent the assembly activity between Xand Z, Zand X, Xthemselves, and Zthemselves.

In reality,


Thus, we determined a general formula for protein assembly. Cytoskeletal protein and protein complexes in the signaling cascade can be described using this formula.

2.7. Linearity of cofactor supply to the assembly system

While protein assembly is a non-linear reaction involving a complicated set of reactions or assembly steps, the supply of cofactor is simply given by the linear kinetic rate items as shown in Eqs. (36, 37, 40 and 41). This means that the supply rate will be essentially be an order parameter of the assembly system and is controllable by altering the concentration of the cofactor (Figure 2). Accordingly, parameter pis variable in the numerical simulation, as described below.

Figure 2.

Scheme of monomer interaction. Individual globules represent monomersX•, and Z ○, released species. The supply of the cofactor is kept constant and inactive cofactor is released continuously. The differential coefficientsa,b,c, anddindicate the interaction activity betweenXandZ.The differential coefficients are given inEq. (39). The interaction activity betweenXis higher and therefore the diffusion rate ofXbetweenXbecomes slower; the interaction activity betweenXandZis lower, and therefore the diffusion rateof X/ZbetweenZ/Xbecomes slightly slower; the interaction activity between Z andZis negligible and therefore the change in the diffusion rate ofZthroughZis negligible.

3. Calculus simulation of concentration oscillation

3.1. Oscillation of monomer concentration fluctuation

In actual simulation of protein assembly, numerical calculation was performed over a sufficiently long period to evaluate the trend in system behavior.

Simulation: A simulation was performed using Mathematica® version 8 (Wolfram Research, Champaign, IL, USA).

Simulation was performed with the notation in Eqs. (42) and (43), in which p is (a) 0.8, (b) 0.81, and (c) 1.00 (Figure 3).

Below is the simulation program cord using Mathematica ver. 8 when p = 0.8:

D1 = 0.27, k2 = 0.00035, a = 790, b = 650, c = 105, d = 105, e = 105, f = 105, p = 0.8, Dxx = 155, Dxz = 155, W = 1.

Figure 3.

Numerical simulation of protein assembly. Diffusion of active cofactor binding signaling molecule (X) and inactive cofactor binding signaling molecule (Z). The upper graph shows two parametric plots of X and Z. Red and blue lines in the lower graph represent the concentrations of X and Z, respectively. The horizontal axis represents time (s) (0 ≤ t ≤ 200) and vertical axis represents the concentrations of X and Z, respectively. When p exceeds 0.80, chaos-like oscillation is observed. Mathematica version 8 was used.

X = k2/D1.

Z = (k2 (D1^2 W + Dxx k2))/(D1 (D1 p - Dxz k2)).

NDSolve[{x’[t] == − (W (D1 - a X) + 2 X Dxx + Dxz Z) x[t] + (W a - Dxx + 2 c X + e Z) x[t]^2 + (p - Dxz X - b X - d X^2 - f X Z) z[t] - (Dxz + W b - e X + f Z) x[t] z[t] - (f X) z[t]^2,

z’[t] == (2 X Dxx + Dxz Z - c X^2 - e X Z) x[t] + (Dxx - 2 c X - e Z) x[t]^2 + (Dxz + 2 X d - e X + f Z) x[t] z[t] + (Dxz X - p + d X^2 + f X Z) z[t], x[0] == 0.000001, z[0] == 0.000001}, {x, z}, {t, 0, 30,000}, MaxSteps - > 50,000].

g001 = Plot[{X + x[t]} /. %, {t, 0, 200}, PlotRange - > All, PlotStyle - > {RGBColor[1, 0, 1]}, PlotRange - > ALL].

g002 = Plot[{Z + z[t]} /. %%, {t, 0, 200}, PlotRange - > All, PlotStyle - > {RGBColor[0, 1, 1]}, PlotRange - > All].

g003 = ParametricPlot[Evaluate[{X + x[t], Z + z[t]} /. %%%], {t, 0, 2000}, PlotRange - > All, AxesLabel - > {“X”, “Z”}] Show[g001, g003, AxesLabel - > {“t”, “X, Z”}].

As a result, regular oscillation with fluctuations in the amplitude and frequency can be illustrated in the plot following the above calculation. Plots on the right side that the oscillation becomes definite within a limit-cycle when p = 1.5.

3.2. Oscillation frequency

In a previous study [6], an interesting relationship was observed between the average frequency of simulated oscillation of the monomeric proteins and difference p − pc. The frequency was nearly equivalent, but with irregular fluctuations, except for during the initial phase. The relationship between average frequency < f > and ppcis given by:


These formulae imply that the amplitude of monomeric protein fluctuation provides information regarding the outside alteration of the cofactor. Thus, outside alteration is transformed inside into the information of assembly.

4. Center manifold theory (CMT)

4.1. Center manifold formulae

For stability analysis of the nonlinear dynamics in protein assembly, the center manifold theory (CMT) for non-linear dynamic biological systems has been applied. Simulation is oriented to analyze the behavior around critical values of the order parameter. CMT has been applied to the Lotka-Volterra model of the predator-prey system to provide important simulation results [33, 34]. However, the CMT can be applied to the protein assembly model. For stability analysis around the critical point, Eqs. (40) and (41) were formulated. When pis equivalent to pc, the Jacobian matrix Lfor (x, z) is given using the linear coefficients of (x, z) in Eqs. (40) and (41):


Particularly, the function f(p) represents the input and output of the cofactor that is the order parameter. Using the eigenvectors of L, (l1, l2), coordinate transformation into uand vis performed as follows:


The above formulae are subsequently set as


The center manifold around the critical point (p = pc) is then given as follows:


The effect of changing pand ε(p = pc) is analyzed using the center manifold around the critical point of the system. Subsequently


Using Eqs. (49) and (50)


Solving Eq. (52) gives the coefficients of aiin Eq. (50): a3 = a7 = 0. Substituting uin Eq. (51) given by νand εinto fv(u, v) in Eq. (47), the kinetic stability equation is given for fluctuation νusing the coefficients ni(i = 1, …, 7) as follows:


Independently of the numerical values in Eq. (53)


Using this result, we have


By setting the left-hand side of Eq. (55) equivalent to zero, a Hopf-bifurcation of the given system is shown


Further approximate solutions to Eq. (56) are given as


Thereafter, the solution increases as the concentration of ATP/GTP increases. From Eq. (51), uis formulated using an arbitrary coefficient c


This result implies that the fluctuation has two different values for the amplitude of an oscillation.

5. Prospects of protein analysis assembly

Our aim in this review was to evaluate the association between the non-linearity in protein diffusion and assembly. Here, we reviewed studies of protein assembly or interactions [1, 2, 3, 4] and performed mathematical analysis of the model in addition to numerical simulations. Because the assumptions of the model are minimal, the simulation provides insight into assembly. The results are summarized as follows: (i) the non-linear kinetic equations including only two independent parameters may reveal dynamic behavior in the fluctuation of the monomer concentration, (ii) the increase near the critical concentration of the cofactor induces oscillations in amplitude and frequency; and (iii) center manifold analysis predicts the stability of the model system near the critical concentration, showing bifurcation with respect to the cofactor supply value. The behavior of the system shown in the simulation indicates that the concentration change information of a cofactor outside the system is transduced into another type of information, e.g., frequency of the concentration oscillation of the monomer. A small increase in the outside cofactor concentration induces an oscillation change inside the monomeric protein, which may be crucial for responding to transformations in the outside environment. Such a trajectory in the observed oscillation resembles a limit-cycle-like in the well-known two-parametric Lorenz model [13, 14].

Previous systems biology models did not focus on the diffusion process of a protein in the cytoplasm. Non-linearity in the process is critical and essential to protein assembly. Before considering a set of simultaneous kinetic equations, non-linearity in the diffusion process should be considered, as de novonucleation is negligible compared with the reaction of the monomer and oligomer. The interaction between assembly-active monomer proteins attenuates the diffusion rate in a non-linear fashion because they can assemble, which inevitably yields non-linearity. As shown in this review, CMT is useful for reducing the parameters of detailed stability analysis around the critical state.


6. Conclusion

Protein interactions play an important role in various biological activities at the cell level. Although protein diffusion is a rate-limiting-step, cell behavior is orchestrated by protein interactions, and signaling transduction, the cytoskeleton, and cell motility are dynamically altered in processes similar to “phase transitions” in inorganic chemical reactions. This review provides a model of phase transition affected by minimal changes in cofactor concentration; but oscillation and bifurcation are inducible by the simple model. Systems biology multi-parametric analysis remains important; however, a simple model sufficiently can better illustrate oscillations in protein concentration with a limit-cycle. Outside alteration such as cofactor concentration change is transformed inside into the information of assembly. Stability analysis using the CMT is a simple method for understanding protein-interacting systems.


This work was supported by the Ministry of Education, Culture, Sport, Science and Technology, Japan, under the project name “Synergy of Fluctuation and Structure.” Project number 2502. The funders had no role in the study design, data collection, and analysis; decision to publish; or preparation of the manuscript.

Conflicts of interest

The authors have declared no conflicts of interest.

© 2017 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Tatsuaki Tsuruyama (December 20th 2017). Non-Linear Kinetic Analysis of Protein Assembly Based on Center Manifold Theory, Kinetic Theory, George Z. Kyzas and Athanasios C. Mitropoulos, IntechOpen, DOI: 10.5772/intechopen.70750. Available from:

chapter statistics

708total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Plasma Kinetic Theory

By Kashif Chaudhary, Auwal Mustapha Imam, Syed Zuhaib Haider Rizvi and Jalil Ali

Related Book

First chapter

Introductory Chapter: Nanomaterials in the 2020s

By George Z. Kyzas and Athanasios C. Mitropoulos

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us