Open access peer-reviewed chapter

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

Written By

Tatsuaki Tsuruyama

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

DOI: 10.5772/intechopen.70750

From the Edited Volume

Kinetic Theory

Edited by George Z. Kyzas and Athanasios C. Mitropoulos

Chapter metrics overview

1,013 Chapter Downloads

View Full Metrics


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 instability in tubulin polymerization has been extensively investigated. Dynamic instability signifies 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, X can reversibly associate with the oligomer:

X + W W kinetic rate coefficients k 1 E1

Subsequently, Z is released from the oligomer:

W Z E2

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

Z + P X + I k 3 E3

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

X + X X + Z k 4 E4
X + Z 2 Z k 5 E5

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 P is 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

P V = D 0 exp βcρ V e / 1 cV e E6

where Ve is 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, D depends on the size and shape of the protein and likely other factors

D = D 0 V e p V dV = D 0 exp βcρ V e / 1 cV e E7

At low concentrations of the macromolecule

D = D 0 1 βcρ V e / 1 cV e D 0 1 + 1 β cV e E8

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

Figure 1.

Self-diffusion rate constants D for concentration of protein, c, with theoretical curve from Eq. (8) with β = 3.0. Near the point x ~ 0, the diffusion coefficient obeys D = 1 − 2c. The vertical axis represents the diffusion coefficient D × 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

η / η 0 = exp ρ η / 1 ρ η k / v E9

where η0 is the solution viscosity at infinite dilution, η is the intrinsic viscosity of the solution, c is the protein solution density, and k/v is 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 X towards oligomer W multiplied by the area of the spherical surface of radius R of the reactive end of polymer or oligomer

Rate of reaction = 4 π R 2 J E10

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

J X = D X d X r dr at r = R = D X X R E11

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

v = 4 πRD X X E12

The rate of the diffusion-controlled reaction is equal to the average flow of X molecule to all W molecules. Accordingly, the global flow of all X to W is 4π R*DXNA XW. Similarly, the flow of all W to X is 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 X to W is given by

Addition rate = 4 πR k X o DXW E13

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

Addition rate = 4 πR k X o D X XW = k X D X XW E14


k X = 4 πR k X o E15

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

dc i dt = k j D i c i + f c i E16

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

J = D i dc j dr + 1 k B T c j r dr E17

In the solvent including two chemical species, the relative movement is related to the local concentration gradient of cj and potential energy. When the gradient is described using Eq. (17), monomer X moves across the sphere with radius R surrounding W and can be described by

J = D X + D W dX dr + 1 k B T X r dr E18

when X and the active site on the oligomer or polymer W interact 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

dX dt = kXW = 4 π R 2 J E19

At the steady state, flux across the sphere with a radius, or a shape parameter, r is constant for any values. Accordingly, R in Eq. (19) is replaced with r using Eq. (18)

dX dt = 4 π R 2 D X + D W dX dr + 1 k B T X r dr E20

By integration and rearrangement of Eq. (20)

kX = 4 π D X + D W γ X = X R X d X r exp ϕ r / k B T E21

and here

γ 1 = R exp ϕ r / k B T r 2 dr E22

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

k = 4 π D X + D W γ X X X R exp ϕ r / k B T = 4 π D X + D W γ 1 + 4 π D X + D W γ / k R exp ϕ r / k B T 4 π D X γ 1 + 4 π D X γ / k R exp ϕ r / k B T E23

Here, Xr=R = (k/kR) X. In calculating the above formula, fluctuations of diffusion coefficients are neglected when the concentration of X is 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 X and W; Xr=R is set to zero. Then, k is given by Eq. (23):

k = 4 π D X R E24

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

D = k B T 6 πrη T E25

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

D = k B T η 1 N A v 1 M 1 c 1 1 + 2 A 1 M 1 c 1 + D 0 1 N A v 1 M 1 c 1 1 + 2 A 1 M 1 c 1 + E26

where T is the temperature of the solution, kB is 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 NA is 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):

a ij D i c j = D 0 i 2 A j M j N A v j M j E27

where vj is 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, X and Z

D X X Z = k B T η X 1 N A v X M X X N A v Z M Z Z 1 + 2 A X M X X + 2 A Z M Z Z + D Z X Z = k B T η Z 1 N A v X M X X N A v Z M Z Z 1 + 2 A X M X X + 2 A Z M Z Z + E28

where vX and vZ are the partial specific volumes of X and Z with molecular weights MX and MZ, respectively. AX and AZ are the second virial coefficients.

2.5. Fluctuation of diffusion coefficient

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

X = X e + x E29
Z = Z e + z E30

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

x ̇ + y ̇ + z ̇ = 0 E31

The dependency of the diffusion rate refers to the high interaction activities of X and low interaction activities of Z. An increase in X contributes to a decrease in the diffusion coefficients DX and DZ in fluctuation items because of the higher interaction activity that reduces diffusion; in contrast, an increased Z contributes to increased diffusion coefficients because the interaction between increased Z induces 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:

X ̇ = k 1 D 1 WX + k 3 PZ k 4 D 4 X 2 k 5 D 5 XZ E32
W ̇ = k 1 D 1 WX k 2 D 2 W 0 E33
Z ̇ = k 2 D 2 W k 3 PZ + k 4 D 4 X 2 + k 5 D 5 XZ E34


D 1 D X + D W 2 D X 2 , D 4 D X , D 5 D X + D Z 2 E35

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

X ̇ = k 1 ' WX + k 3 PZ k 4 ' X 2 k 5 ' XZ E36
Z ̇ = k 2 ' W k 3 PZ + k 4 ' X 2 + k 5 ' XZ E37

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

X e = k 2 k 1 , Z e = k 2 k 4 4 k 2 + k 1 2 W k 1 k 1 P k 2 k 5 E38

Here, the fluctuation of diffusion coefficients are given by

a = k 1 ' x , b = k 1 ' z , c = k 4 ' x , d = k 4 ' z , e = k 5 ' x , f = k 5 ' z E39

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

x ̇ = f x x + f z + p z + f xx x 2 + f xz xz + f zz z 2 E40
z ̇ = f x ' x + f z ' f z p z + f xx ' x 2 + f zx ' xz + f zz ' z 2 E41

Accordingly, the overall behavior of the kinetics of protein assembly is given by the monomeric kinetics of x and z. fxz(), fzx(), fxx(), and fzz() represent the assembly activity between X and Z, Z and X, X themselves, and Z themselves.

In reality,

x ̇ = W k 1 aX e + 2 k 4 X e + k 5 Z e x + Wa k 4 + 2 cX e + eZ e x 2 + p bX e k 5 X e dX e 2 fX e Z e z k 5 + Wb eX e + f e Z e xz fX e z 2 E42
z ̇ = 2 k 4 X e + k 5 Z e cX e 2 eX e Z e x + k 4 2 cX e eZ e x 2 + p + k 5 X e + dX e 2 + fX e Z e z + k 5 + 2 dX e eX e + f e Z e xz + fX e z 2 E43

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 p is variable in the numerical simulation, as described below.

Figure 2.

Scheme of monomer interaction. Individual globules represent monomers X•, and Z ○, released species. The supply of the cofactor is kept constant and inactive cofactor is released continuously. The differential coefficients a, b, c, and d indicate the interaction activity between X and Z. The differential coefficients are given in Eq. (39). The interaction activity between X is higher and therefore the diffusion rate of X between X becomes slower; the interaction activity between X and Z is lower, and therefore the diffusion rate of X/Z between Z/X becomes slightly slower; the interaction activity between Z and Z is negligible and therefore the change in the diffusion rate of Z through Z is 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 ppc is given by:

< f > 0.0256 ln p p c + 0.1407 E44

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 p is equivalent to pc, the Jacobian matrix L for (x, z) is given using the linear coefficients of (x, z) in Eqs. (40) and (41):

L = f x f z + f p c f x ' f z ' F p c E45

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 u and v is performed as follows:

d dt x z = d dt l 1 l 2 u v E46
d dt u v = d dt l 1 l 2 1 x z E47

The above formulae are subsequently set as

du / dt = f u u v E48
dv / dt = f v u v E49

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

u = h ε v = a 1 v 2 + a 2 + a 3 ε 2 + a 4 v 3 + a 5 v 2 ε + a 6 v ε 2 + a 7 ε 3 + O ε 4 E50

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

u = dv / dt h u ε / u + / dt h u ε / ε = 2 a 1 v + a 2 ε f u u v E51

Using Eqs. (49) and (50)

2 a 1 v + a 2 ε f u u v = a 1 v 2 + a 2 + a 3 ε 2 + a 4 v 3 + a 5 v 2 ε + a 6 v ε 2 + a 7 ε 3 + O ε 4 E52

Solving Eq. (52) gives the coefficients of ai in Eq. (50): a3 = a7 = 0. Substituting u in 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:

dv / dt = n 1 v 2 + n 2 + n 3 ε 2 + n 4 v 3 + n 5 v 2 ε + n 6 v ε 2 + n 7 ε 3 + O ε 4 E53

Independently of the numerical values in Eq. (53)

n 3 , n 6 , n 7 = 0 E54

Using this result, we have

dv / dt = n 1 v 2 + n 2 + n 4 v 3 + n 5 v 2 ε + O ε 4 E55

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

v = 0 , n 1 n 5 ε ± n 1 + n 5 ε 2 4 n 2 n 4 ε 2 n 4 E56

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

v = 0 , , n 1 / n 4 E57

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

u 0 , c n 1 / n 4 2 E58

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 novo nucleation 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.


  1. 1. Gorecki J, Gorecka JN, Nowakowski B, Ueno H, Tsuruyama T, Yoshikawa K. Sensing parameters of a time dependent inflow with an enzymatic reaction. In: Adamatzky A, editor. Advances in Unconventional Computing, Emergence, Complexity and Computation. Vol. 23; 2017. p. 85-104
  2. 2. Tsuruyama T. A model of cell biological signaling predicts a phase transition of signaling and provides mathematical formulae. PloS One. 2014;9(7):e102911
  3. 3. Tsuruyama T. Kinetic stability analysis of protein assembly on the center manifold around the critical point. BMC Systems Biology. 2017;11(1):13
  4. 4. Ueno H, Tsuruyama T, Nowakowski B, Gorecki J, Yoshikawa K. Discrimination of time-dependent inflow properties with a cooperative dynamical system. Chaos. 2015;25(10):103115
  5. 5. Di Camillo B, Carlon A, Eduati F, Toffolo GM. A rule-based model of insulin signalling pathway. BMC Systems Biology. 2016;10(1):38
  6. 6. Babu CVS, Song EJ, Yoo YS. Modeling and simulation in signal transduction pathways: A systems biology approach. Biochimie. 2006;88(3–4):277-283
  7. 7. Selimkhanov J, Taylor B, Yao J, Pilko A, Albeck J, Hoffmann A, Tsimring L, Wollman R. Systems biology. Accurate information transmission through dynamic biochemical signaling networks. Science. 2014;346(6215):1370-1373
  8. 8. Zhao M, Kong L, Qu H. A systems biology approach to identify intelligence quotient score-related genomic regions, and pathways relevant to potential therapeutic treatments. Scientific Reports. 2014;4:4176
  9. 9. Hazra P, Inoue K, Laan W, Hellingwerf KJ, Terazima M. Tetramer formation kinetics in the signaling state of AppA monitored by time-resolved diffusion. Biophysical Journal. 2006;91(2):654-661
  10. 10. Wu Z, Wang HW, Mu W, Ouyang Z, Nogales E, Xing J. Simulations of tubulin sheet polymers as possible structural intermediates in microtubule assembly. PloS One. 2009;4(10):e7291
  11. 11. VanBuren V, Cassimeris L, Odde DJ. Mechanochemical model of microtubule structure and self-assembly kinetics. Biophysical Journal. 2005;89(5):2911-2926
  12. 12. Hamon L, Savarin P, Curmi PA, Pastré D. Rapid assembly and collective behavior of microtubule bundles in the presence of polyamines. Biophysical Journal. 2011;101(1):205-216
  13. 13. Wattis JAD, Coveney PV. Mesoscopic models of nucleation and growth processes: A challenge to experiment. Physical Chemistry Chemical Physics. 1999;1:2163-2176
  14. 14. Carlsson AE. Model of reduction of actin polymerization forces by ATP hydrolysis. Physical Biology. 2008;5(3):036002
  15. 15. Brooks FJ, Carlsson AE. Actin polymerization overshoots and ATP hydrolysis as assayed by pyrene fluorescence. Biophysical Journal. 2008;95(3):1050-1062
  16. 16. Oosawa F, Kasai M. A theory of linear and helical aggregations of macromolecules. Journal of Molecular Biology. 1962;4:10-21
  17. 17. Michaels TC, Garcia GA, Knowles TP. Asymptotic solutions of the Oosawa model for the length distribution of biofilaments. The Journal of Chemical Physics. 2014;140(19):194906
  18. 18. Chretien D, Jainosi I, Taveau JC, Flyvbjerg H. Microtubule’s conformational cap. Cell Structure and Function. 1999;24(5):299-303
  19. 19. Zilberman M, Sofer M. A mathematical model for predicting controlled release of bioactive agents from composite fiber structures. Journal of Biomedical Materials Research. Part A. 2007;80(3):679-686
  20. 20. Oosawa F, Asakura S. Thermodynamics of the Polymerisation of Proteins; New York: Academic Press, 1975. p. 204
  21. 21. Zapperi S, Mahadevan L. Dynamic instability of a growing adsorbed polymorphic filament. Biophysical Journal. 2011;101(2):267-275
  22. 22. Hammele M, Zimmermann W. Modeling oscillatory microtubule polymerization. Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics. 2003;67(2 Pt 1):021903
  23. 23. Chen BS, CC W. On the calculation of signal transduction ability of signaling transduction pathways in intracellular communication: Systematic approach. Bioinformatics. 2012;28(12):1604-1611
  24. 24. Zumsande M, Gross T. Bifurcations and chaos in the MAPK signaling cascade. Journal of Theoretical Biology. 2010;265(3):481-491
  25. 25. Arnal I, Karsenti E, Hyman AA. Structural transitions at microtubule ends correlate with their dynamic properties in Xenopus egg extracts. Journal of Cell Biology. 2000;149(4):767-774
  26. 26. Daga RR, Lee KG, Bratman S, Salas-Pino S, Chang F. Self-organization of microtubule bundles in anucleate fission yeast cells. Nature Cell Biology. 2006;8(10):1108-1113
  27. 27. Kasas S, Cibert C, Kis A, De Los Rios P, Riederer BM, Forro L, Dietler G, Catsicas S. Oscillation modes of microtubules. Biology of the Cell. 2004;96(9):697-700
  28. 28. Bauer KC, Göbel M, Schwab ML, Schermeyer MT, Hubbuch J. Concentration-dependent changes in apparent diffusion coefficients as indicator for colloidal stability of protein solutions. International Journal of Pharmaceutics. 2016;511(1):276-287
  29. 29. Gallagher WH, Woodward CK. The concentration dependence of the diffusion coefficient for bovine pancreatic trypsin inhibitor: A dynamic light scattering study of a small protein. Biopolymers. 1989;28(11):2001-2024
  30. 30. Giavazzi F, Fornasieri A, Vailati A, Cerbino R. Equilibrium and non-equilibrium concentration fluctuations in a critical binary mixture. The European Physical Journal. E, Soft Matter. 2016;39(10):103
  31. 31. Giavazzi F, Savorana G, Vailati A, Cerbino R. Structure and dynamics of concentration fluctuations in a non-equilibrium dense colloidal suspension. Soft Matter. 2016;12(31):6588-6600
  32. 32. Chang X, Wei J. Stability and Hopf bifurcation in a diffusive predator-prey system incorporating a prey refuge. Mathematical Biosciences and Engineering. 2013;10(4):979-996
  33. 33. Zhang X, Zhao H. Bifurcation and optimal harvesting of a diffusive predator-prey system with delays and interval biological parameters. Journal of Theoretical Biology. 2014;363:390-403
  34. 34. Xiao M, Zheng WX, Cao J. Hopf bifurcation of an (n + 1)-neuron bidirectional associative memory neural network model with delays. IEEE Transactions on Neural Networks and Learning Systems. 2013;24(1):118-132

Written By

Tatsuaki Tsuruyama

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