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

Abstract

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.

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\to W\left(\text{kinetic rate coefficients}\to k_1\right)$ E1

Subsequently, Z is released from the oligomer:

$W\to 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\to X+I\left(k_3\right)$ E3

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

$X+X\to X+Z\left(k_4\right)$ E4

$X+Z\to 2Z\left(k_5\right)$ 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 D₀ 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 = D_oP.

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

$P\left(V\right)=D_0\text{exp}\left[-\mathit{\beta c\rho }{V}_e/\left(1-{\mathit{cV}}_e\right)\right]$ E6

where V_e 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{\int }_{V_e}^{\infty }p\left(V\right)\mathit{dV}=D_0\text{exp}\left[-\mathit{\beta c\rho }{V}_e/\left(1-{\mathit{cV}}_e\right)\right]$ E7

At low concentrations of the macromolecule

$D=D_0\left(1-\mathit{\beta c\rho }{V}_e/\left(1-{\mathit{cV}}_e\right)\right)\simeq D_0\left(1+\left(1-\beta \right){\mathit{cV}}_e\right)$ 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 × 10⁷ (cm²/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

$\eta /{\eta }_0=\text{exp}\left(\rho \left(\eta \right)/\left(1-\rho \left(\eta \right)\left(k/v\right)\right)\right)$ E9

where η₀ 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

$\text{Rate of reaction}=4\pi 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{\left(\frac{d{\left[X\right]}_r}{\mathit{dr}}\right)}_{\mathit{at}\mathrm{r}=\mathrm{R}}=D_X\frac{\left[X\right]}{R}$ E11

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

$v=4\mathit{\pi 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*D_XN_A XW. Similarly, the flow of all W to X is 4π RD_WN_A WX. Further, using the sum of the diffusion coefficients of the two species, the diffusion coefficient is rewritten as D = (D_X + D_W)/2. Then, the addition rate of X to W is given by

$\text{Addition rate}=4\mathit{\pi R}{k_X}^o\mathit{DXW}$ E13

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

$\text{Addition rate}=4\mathit{\pi R}{k_X}^o{D}_X\mathit{XW}=k_X{D}_X\mathit{XW}$ E14

and

$k_X=4\mathit{\pi R}{k_X}^o$ E15

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

$\frac{{\mathit{dc}}_i}{\mathit{dt}}=k_j{D}_i{c}_i+f\left(c_i\right)$ E16

where f(c_j) represents the sum of the reaction kinetic items of jth species, c_j, 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\left[\frac{{\mathit{dc}}_j}{\mathit{dr}}+\frac{1}{k_{B}T}{c}_j\frac{\mathit{d\varphi }\left(r\right)}{\mathit{dr}}\right]$ E17

In the solvent including two chemical species, the relative movement is related to the local concentration gradient of c_j 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=-\left(D_X+D_W\right)\left[\frac{\mathit{dX}}{\mathit{dr}}+\frac{1}{k_{B}T}X\frac{\mathit{d\varphi }\left(r\right)}{\mathit{dr}}\right]$ 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

$\frac{\mathit{dX}}{\mathit{dt}}=-\mathit{kXW}=-4\pi 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)

$\frac{\mathit{dX}}{\mathit{dt}}=-4\pi R^2\left(D_X+D_W\right)\left[\frac{\mathit{dX}}{\mathit{dr}}+\frac{1}{k_{B}T}X\frac{\mathit{d\varphi }\left(r\right)}{\mathit{dr}}\right]$ E20

By integration and rearrangement of Eq. (20)

$\mathit{kX}=-4\pi \left(D_X+D_W\right)\gamma {\int }_{X=X_R}^{X}d\left\{X_r\text{exp}\left(\varphi \left(r\right)/k_{B}T\right)\right\}$ E21

and here

${\gamma }^{-1}={\int }_R^{\infty }\frac{\text{exp}\left(\varphi \left(r\right)/k_{B}T\right)}{r^2}\mathit{dr}$ E22

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

$\begin{array}{l}k=\frac{4\pi \left(D_X+D_W\right)\gamma }{X}\left\{X-X_R\text{exp}\left(\varphi \left(r\right)/k_{B}T\right)\right\}\\ =\frac{4\pi \left(D_X+D_W\right)\gamma }{1+\left[4\pi \left(D_X+D_W\right)\gamma /k_R\text{exp}\left(-\varphi \left(r\right)/k_{B}T\right)\right]}\\ \approx \frac{4\pi D_X\gamma }{1+\left[4\pi D_X\gamma /k_R\text{exp}\left(-\varphi \left(r\right)/k_{B}T\right)\right]}\end{array}$ E23

Here, X_r=R = (k/k_R) 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; X_r=R is set to zero. Then, k is given by Eq. (23):

$k=4\pi 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=\frac{k_{B}T}{6\mathit{\pi r\eta }\left(T\right)}$ E25

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

$D=\frac{k_{B}T}{\eta }\left(1-\frac{N_A{v}_1}{M_1}{c}_1\right)\left(1+2A_1{M}_1{c}_1+\cdots \right)\triangleq D_0\left(1-\frac{N_A{v}_1}{M_1}{c}_1\right)\left(1+2A_1{M}_1{c}_1+\cdots \right)$ E26

where T is the temperature of the solution, k_B is the Boltzmann constant, and η₁ is the frictional coefficient of the macromolecule in solution. A₁ is the second virial coefficient, v₁ is the partial specific volume of protein with molecular weight M₁, and N_A is Avogadro’s number. The small letter c₁ denotes the concentration of the solute macromolecule. Then, dependency of the diffusion coefficient on the ith component, D_i, is as follows from (26):

$a_{\mathit{ij}}\equiv \frac{\partial D_i}{\partial c_j}=D_{0i}\left(2A_j{M}_j-\frac{N_A{v}_j}{M_j}\right)$ E27

where v_j is the partial specific volume of the polymer with molecular weight M_j.

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

$\begin{array}{l}{D}_X\left(X,Z\right)=\frac{k_{B}T}{{\eta }_X}\left(1-\frac{N_A{v}_X}{M_X}X-\frac{N_A{v}_Z}{M_Z}Z\right)\left(1+2A_X{M}_{X}X+2A_Z{M}_{Z}Z+\cdots \right)\\ D_Z\left(X,Z\right)=\frac{k_{B}T}{{\eta }_Z}\left(1-\frac{N_A{v}_X}{M_X}X-\frac{N_A{v}_Z}{M_Z}Z\right)\left(1+2A_X{M}_{X}X+2A_Z{M}_{Z}Z+\cdots \right)\end{array}$ E28

where v_X and v_Z are the partial specific volumes of X and Z with molecular weights M_X and M_Z, respectively. A_X and A_Z 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

$\dot{x}+\dot{y}+\dot{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 D_X and D_Z 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:

$\dot{X}=-k_1{D}_1\mathit{WX}+k_3\mathit{PZ}-k_4{D}_4{X}^2-k_5{D}_5\mathit{XZ}$ E32

$\dot{W}=k_1{D}_1\mathit{WX}-k_2{D}_{2}W\approx 0$ E33

$\dot{Z}=k_2{D}_{2}W-k_3\mathit{PZ}+k_4{D}_4{X}^2+k_5{D}_5\mathit{XZ}$ E34

Here,

$D_1\triangleq \frac{D_X+D_W}{2}\sim \frac{D_X}{2},D_4\triangleq D_X,D_5\triangleq \frac{D_X+D_Z}{2}$ E35

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

$\dot{X}=-{k_1}^{\text{'}}\mathit{WX}+k_3\mathit{PZ}-{k_4}^{\text{'}}{X}^2-{k_5}^{\text{'}}\mathit{XZ}$ E36

$\dot{Z}={k_2}^{\text{'}}W-k_3\mathit{PZ}+{k_4}^{\text{'}}{X}^2+{k_5}^{\text{'}}\mathit{XZ}$ E37

Here, k₁D₁W = k₁’, k₄D₄W = k₄’, k₅D₅ = k₅’, and p = k₃ P. At the steady state, setting the right hands of Eqs. (32)-(34) equal to zero

$X_e=\frac{{k_2}^{\prime }}{{k_1}^{\prime }},Z_e=\frac{{k_2}^{\prime }\left({k_4}_4^{\prime }{k_2}^{\prime }+{k_1}^{2\prime }W\right)}{{k_1}^{\prime }\left({k_1}^{\prime }P-{k_2}^{\prime }{k_5}^{\prime }\right)}$ E38

Here, the fluctuation of diffusion coefficients are given by

$a=\frac{\partial {k_1}^{\text{'}}}{\partial x},b=\frac{\partial {k_1}^{\text{'}}}{\partial z},c=\frac{\partial {k_4}^{\text{'}}}{\partial x},d=\frac{\partial {k_4}^{\text{'}}}{\partial z},e=\frac{\partial {k_5}^{\text{'}}}{\partial x},f=\frac{\partial {k_5}^{\text{'}}}{\partial z}$ E39

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

$\dot{x}=f_{x}x+\left(f_z+p\right)z+f_{\mathit{xx}}{x}^2+f_{\mathit{xz}}\mathit{xz}+f_{\mathit{zz}}{z}^2$ E40

$\dot{z}={f_x}^{\text{'}}x+\left({f_z}^{\text{'}}-f_z\left(p\right)\right)z+{f_{\mathit{xx}}}^{\text{'}}{x}^2+{f_{\mathit{zx}}}^{\text{'}}\mathit{xz}+{f_{\mathit{zz}}}^{\text{'}}{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,

$\begin{array}{l}\dot{x}=-\left\{W\left({k_1}^{\prime }-{\mathit{aX}}_e\right)+2{k_4}^{\prime }{X}_e+{k_5}^{\prime }{Z}_e\right\}x+\left(\mathit{Wa}-{k_4}^{\prime }+2{\mathit{cX}}_e+{\mathit{eZ}}_e\right)x^2\\ +\left(p-{\mathit{bX}}_e-{k_5}^{\prime }{X}_e-{{\mathit{dX}}_e}^2-{\mathit{fX}}_e{Z}_e\right)z-\left({k_5}^{\prime }+\mathit{Wb}-{\mathit{eX}}_e+f_e{Z}_e\right)\mathit{xz}-{\mathit{fX}}_e{z}^2\end{array}$ E42

$\begin{array}{l}\dot{z}=\left(2{k_4}^{\prime }{X}_e+{k_5}^{\prime }{Z}_e-{{\mathit{cX}}_e}^2-{\mathit{eX}}_e{Z}_e\right)x+\left({k_4}^{\prime }-2{\mathit{cX}}_e-{\mathit{eZ}}_e\right)x^2\\ +\left(-p+{k_5}^{\prime }{X}_e+{{\mathit{dX}}_e}^2+{\mathit{fX}}_e{Z}_e\right)z+\left({k_5}^{\prime }+2{\mathit{dX}}_e-{\mathit{eX}}_e+f_e{Z}_e\right)\mathit{xz}+{\mathit{fX}}_e{z}^2\end{array}$ 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 − p_c. The frequency was nearly equivalent, but with irregular fluctuations, except for during the initial phase. The relationship between average frequency < f > and p − p_c is given by:

$<f>\simeq 0.0256\mathit{ln}\left(p-p_c\right)+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 p_c, the Jacobian matrix L for (x, z) is given using the linear coefficients of (x, z) in Eqs. (40) and (41):

$L=\left(\begin{array}{ll}{f}_x\hfill & f_z+f\left(p_c\right)\hfill \\ {f_x}^{\text{'}}\hfill & {f_z}^{\text{'}}-F\left(p_c\right)\hfill \end{array}\right)$ E45

Particularly, the function f(p) represents the input and output of the cofactor that is the order parameter. Using the eigenvectors of L, (l₁, l₂), coordinate transformation into u and v is performed as follows:

$\frac{d}{\mathit{dt}}\left[\begin{array}{l}x\hfill \\ z\hfill \end{array}\right]=\frac{d}{\mathit{dt}}\left[{\mathbf{l}}_{\mathbf{1}}{\mathbf{l}}_{\mathbf{2}}\right]\left[\begin{array}{l}u\hfill \\ v\hfill \end{array}\right]$ E46

$\frac{d}{\mathit{dt}}\left[\begin{array}{l}u\hfill \\ v\hfill \end{array}\right]=\frac{d}{\mathit{dt}}{\left[\begin{array}{ll}{\mathbf{l}}_{\mathbf{1}}\hfill & {\mathbf{l}}_{\mathbf{2}}\hfill \end{array}\right]}^{-1}\left[\begin{array}{l}x\hfill \\ z\hfill \end{array}\right]$ E47

The above formulae are subsequently set as

$\mathit{du}/\mathit{dt}=f_u\left(u,v\right)$ E48

$\mathit{dv}/\mathit{dt}=f_v\left(u,v\right)$ E49

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

$u=h\left(\epsilon ,v\right)=a_1{v}^2+a_2\mathit{v\epsilon }+a_3{\epsilon }^2+a_4{v}^3+a_5{v}^2\epsilon +a_{6}v{\epsilon }^2+a_7{\epsilon }^3+O\left({\epsilon }^4\right)$ E50

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

$u=\left(\mathit{dv}/\mathit{dt}\right)\partial h\left(u,\epsilon \right)/\partial u+\left(\mathit{d\epsilon }/\mathit{dt}\right)\partial h\left(u,\epsilon \right)/\partial \epsilon =\left(2a_{1}v+a_2\epsilon \right)f_u\left(u,v\right)$ E51

Using Eqs. (49) and (50)

$\left(2a_{1}v+a_2\epsilon \right)f_u\left(u,v\right)=a_1{v}^2+a_2\mathit{v\epsilon }+a_3{\epsilon }^2+a_4{v}^3+a_5{v}^2\epsilon +a_{6}v{\epsilon }^2+a_7{\epsilon }^3+O\left({\epsilon }^4\right)$ E52

Solving Eq. (52) gives the coefficients of a_i in Eq. (50): a₃ = a₇ = 0. Substituting u in Eq. (51) given by ν and ε into f_v(u, v) in Eq. (47), the kinetic stability equation is given for fluctuation ν using the coefficients n_i (i = 1, …, 7) as follows:

$\mathit{dv}/\mathit{dt}=n_1{v}^2+n_2\mathit{v\epsilon }+n_3{\epsilon }^2+n_4{v}^3+n_5{v}^2\epsilon +n_{6}v{\epsilon }^2+n_7{\epsilon }^3+O\left({\epsilon }^4\right)$ E53

Independently of the numerical values in Eq. (53)

$n_3,n_6,n_7=0$ E54

Using this result, we have

$\mathit{dv}/\mathit{dt}=n_1{v}^2+n_2\mathit{v\epsilon }+n_4{v}^3+n_5{v}^2\epsilon +O\left({\epsilon }^4\right)$ E55

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

$v=0,\frac{-n_1-n_5\epsilon \pm \sqrt{{\left(n_1+n_5\epsilon \right)}^2-4n_2{n}_4\epsilon }}{2n_4}$ E56

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

$v=0,\approx \mathit{c\epsilon },-n_1/n_4$ E57

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

$u\approx 0,c{\left(n_1/n_4\right)}^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.

Acknowledgments

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.