Mathematical modeling and analysis of the pathway network consisting of symmetrical complexes with N monomers, like the activation of MMP2
MMathematical modeling and analysis of the pathway networkconsisting of symmetrical complexes with N monomers, like theactivation of MMP2.
Keiko Itano ∗ July 5, 2019
Abstract
The activation of matrix metalloproteinase 2 (MMP2) is a crucial event during tumor metastasisand invasion, and this pathway network consists of 3 monomers. The pathway network of the ac-tivation obeys to a set of specified reaction rules. According to the rules, the individual moleculeslocalize in a particular order and symmetrically around a homodimer following the formation of thatdimer.We generalized the homodimer pathway network obeying to similar reaction rules, by changingthe number of monomers involved in this pathway from 3 to N. At the previous work, we found themolecules in the pathway network are classified to some reaction groups. We derived the law of massconservation between the groups. Each group concentration converges to its equilibrium solution.Using these results, we derive the concentrations of the complexes theoretically and reveal that eachcomplex concentration converges to its equilibrium value. We can say the pathway network withhomodimer symmetric form complexes is asymptotic stable and identify the regulator parameter of thetarget complex in the network. Our mathematical approach may help us understand the mechanismof this type pathway network by knowing the background mathematical laws which govern this typepathway network.
Matrix metalloproteinases (MMPs) represent a family of endopeptidases that are responsible forthe degradation of many components of extracellular matrix (ECM). The process of ECM degradationplays an important role during tumor metastasis and invasion, and this pathway network consists of3 monomers: matrix metalloproteinase 2 (MMP2), tissue inhibitor of metalloproteinase(TIMP2), andmembrane type-1 matrix metalloproteinase(MT1-MMP). Sato et al. (1994) demonstrated that MMP2activation occurs upon the formation of (MMP2- TIMP2 - MT1-MMP - MT1-MMP) complex [1] (Fig.1).MMP inhibitors for cancer treatments have been developed, but their many side effects have hinderedthe use of these inhibitors. Therefore, the suppression of this complex is necessary in order to developimproved treatments, and the mathematical modeling mechanisms involved in the activation of MMP2may represent a very useful approach in anti-cancer drug development. Several mathematical approachesof the activation pathway of MMP 2 have been proposed. Hoshino et al. (2012) proposed a computationalmodel of this (MMP2 - TIMP2- MT1-MMP - MT1-MMP) complex, taking into consideration the transientdynamics of the MMP2 activation [3]. Saito et al. (2012) investigated this computational model, and, usingthe computational simulations, found a novel drug target useful for the suppression of (MMP2-TIMP2- ∗ Osaka University, [email protected] a r X i v : . [ q - b i o . B M ] M a r odel Figure 1. MT1-MMP complex activates MMP2. When the (MMP2 - TIMP2 - MT1-MMP - MT1-MMP) complex is formed,the interactions between MMP2 and TIMP2 are disrupted, which leads to the activation of MMP2.
MT1-MMP - MT1-MMP) complex during the MMP2 activation process [4].The ordinary differential equations (ODEs), as the molecular network evolution equations, have toomany terms to be solved theoretically. Kawasaki et al. showed that the MMP2 activation pathwaynetwork can be divided into several groups, and that ODEs of the MMP2 activation network are solvableand all complexes have explicit solutions[5]. Itano and Suzuki extended the 3-monomer model to the Nmonomer model and showed that the monomers in the extended pathway network can be classified intoreaction groups. By grouping, the N ( N +1) molecule ODEs are summed into the N group ODEs accordingto the law of mass conservation. The group ODEs are solved explicitely and have their equilibria. [6]The aim of this study was to understand the mechanism of a biochemical reaction by generalizing apathway network and solving the ODEs, in order to elucidate the relationships between the moleculesinvolved in this pathway, to investigate the pathway network regulation, and to identify theoreticallyparameters important for the regulation of this network. The generalization of the pathway network,including N monomers, allows a better understanding of the biochemical reaction mechanism.In the complexes formed during MMP2 activation, the monomers localize in a particular order, sym-metrically around the MT1-MMP homodimer. We found that the formation of the molecular complexesin the networks depends on the chemical reaction rules of the network and decides the pathway networkevolution.According to their role in the complex formation, the molecules of the network are classified into severalreaction groups, with N mass conservation laws between the groups. Using these mass conservation laws,the N(N+1) molecule ODEs can be simplified to N group ODEs, which have strict solutions and theequilibrium solutions. The equilibrium solutions of the molecule ODEs are calculated theoretically withthe group solutions as their own upper values.In this paper, we introduce MMP2 activation model as an example of a pathway network that showsa symmetrical complex formation, ( b N − · · · − b − b − · · · − b N ) type complex formation. Most polymersform circular symmetric complexes, with the b N − b N connections, because this ring formation allows for amore stable formation of the complexes. However, ( b N −· · ·− b − b −· · ·− b N ) formation type symmetricpolymers exists, and it is possible to apply the results obtained in this study to the pathway networksin which molecules have similar symmetric complex formation, such as the formation of histone octamer[9], autophagy complex (Atg8-Atg3-Atg12-Atg5-Atg16-Atg16-Atg5-Atg12)[11] [12] [13], semaphorin andplexin complex [10], Kai3-Kai2-Kai1 complex, immunoglobulin complex, and others (Fig. 2)Shinar and Feinberg investigated the robustness of network depending on its structure. They divided c) (a) (b) (d) Figure 2. Examples of the symmetric formation of biological ( b N − · · · − b − b − · · · − b N ) type complexes. (a) Histoneoctamer wrapped with DNA takes a (H2B-H2A-H4-H3-H3-H4-H2A-H2B) form. This complex plays an important roleduring transcription. (b) Semaphorin and plexin complex takes a (PlxnA2-Sema-Sema-PlxnA2) form. (c) (Atg12-Atg5-Atg16-Atg16-Atg5-Atg12) complex reacts with Atg3-Atg8 complex and forms (Atg8-Atg3-Atg12-Atg5-Atg16-Atg16-Atg5-Atg12). This complex formation is similar to that of (MMP2-TIMP2-MT1-MMP-MT1-MMP). a pathway network into different modules, and showed the robustness of the network depends on therobustness of all reaction modules [7]. Mochizuki and Fiedler describe the sensitivity of the flux responseto a change in the reaction rate. They investigated the network structure in order to obtain the responsesensitivity. For a pathway network with a layered structure, they divided the network to motifs andobtained a local flux response. Afterward, they obtained the global flux response sensitivity [8]. Inthe network structure approach, a network is divided into several modules, and the local modules areindependent of each other. The local module parameters or characteristic are obtained independently,while the total network characteristics are investigated globally.Our approach is similar to the described approaches from the perspective that we classify the pathwaynetwork to some modules and clarify the structure. The extended N monomer pathway network has alayered structure because the complexes are classified into several groups according to their formation.We solved all molecular concentrations, using group solutions obtained previously [6], where we solvedgroup solutions locally. Here, we described and solved the complex solutions globally. We generalized a 3-monomer pathway network, obtaining the N monomer pathway network. Thereaction rules are similar to those in the 3-monomer MMP2 activation network. b monomers react eachother and make homodimer. Other monomers connects with monomers in a specified order. We showthe reaction rules as follows. b monomer reacts with b monomer and makes a homodimer b − b . b + b → b b , ( k : reaction rate constant and l : the dissolution rate constant ) • b i monomer reacts with b i − or b i +1 monomer. ( i = 2 , · · · , N − b i +1 + b i → b i +1 b i , ( k i +1 :the reaction rate constant , and l i +1 :the dissolution rate constant)According to the reaction rules, the pathway network has N ( N + 1) types of complexes as follows. b N , b N b N − , . . . , b N b N − · · · b b , b N · · · b b , b N · · · b b b , . . . , b N · · · b b · · · b N b N − , b N − b N − , . . . , b N − · · · − b , b N − · · · b b , . . . , b N − · · · b b · · · b N − ... b , b b , b b b ,b , b b N monomer pathway network and its structural features The generalization helps us understand the structural features dependent on to the complex formations.We determined that different groups of complexes react with each other, in order to produce a new typeof complex, and this grouping depends on the formation of complex and the type of the edge monomer.
The complexes belong to the reaction groups, and we determined mass conservation relationships betweenthe reaction groups. N ( N + 1) ODEs for complex evolution analyses were aggregated to the N ODEsfor reaction groups. The N reaction group ODEs are independent of each other and quadratic with onevariable, and therefore, solvable. We showed that reaction group ODEs have strict solutions and thesolutions converge to the equilibrium state. The reaction group solutions were obtained. First, we obtained a priori upper bound of the concentrationof each complex in the group, and afterward, we show that each complex ODE is written as follows: dX l,m ( t ) dt = − A l,m ( t ) X l,m ( t ) + f l,m ( t ) , (1)where variables A l,m ( t ) and f l,m ( t ) can be written with reaction group variable ξ m ( t ) and reaction rateconstants k m and l m ( m =1, · · · , N). Reaction group variable ξ m ( t ) is strictly derived. The solution of thisequation is: X l,m ( t ) = e − (cid:82) t A l,m ( s ) ds + (cid:90) t e − (cid:82) s A l,m ( u ) du f l,m ( s ) ds. (2)This shows that each complex concentration has a solution, and the N ( N + 1) complex ODEs areintegrable. Using the equation (2), we obtained the equilibrium solution. Therefore, the equilibriumsolution of X ∗ l,m = lim t →∞ X l,m ( t ) is: X ∗ l,m = f ∗ l,m A ∗ l,m . (3) B ・・・ B N-2 B N-1 B N
1 b b ・・・ ・・・ b N-1 b N
2 b b b b ・・・ ・・・ b N-1 b N-2 b N b N-1 ・・・ b b b ・・・ ・・・ b N-1 b N-2 b N-3 ・・・ ・・・ b b b b ・・・ ・・・ ・・・ ・・・ N-1 ・・・ ・・・ b N-1 b N-2 ・・・ b b N b N-1 ・・・ b N ・・・ ・・・ b N-1 b N-2 ・・・ b b b N b N-1 ・・・ b b N+1 ・・・ ・・・ b N b N-1 ・・・ b b b ・・・ ・・・ ・・・ ・・・ ・・・ b N-1 ・・・ b b ・・ b N-2 ・・・ b N-1 ・・・ b b ・・・ b N-2 b N-1 b N b N-1 ・・・ b b ・・ b N-2 b N b N-1 ・・・ b b ・・・ b N-1 b N b N-1 ・・・ b b ・・・ b N-1 b N 複合体の形状がネット ワーク構造・数理モデル に大きく影響 b b を基点とする Figure 3. Diagram molecules belonging to the N monomer network. Rows represent the number of monomers in thecomplex. The column is the biggest monomer index number. The molecules are expressed with the complex formationconsisting N monomers, like ( b m − · · · − b − b − · · · − b l − m ). N monomer network Our aim was to understand the behavior of the investigated network.The N monomer network has N ( N +1) types of molecules. (Fig.3-Fig.4). The molecular concentrationis expressed as X l,m ( t ), where m is the biggest index number of a monomer in the molecule, while l isthe length of the molecule. We classified the molecules to several reaction groups according to the resulting reactions of the molecule X l,m . (Fig.5)First, we classified the molecules to the N groups according to the biggest monomer index numberin the molecule. The monomer with the biggest index number is the furthest from the central monomer b . We defined a reaction group as B i ( i = 1 , ..., N ). Group B i consists of molecules, X ,i : ( b i ), X ,i :( b i − b i − ), ..., X i,i : ( b i −· · ·− b ), X i +1 ,i : ( b i −· · ·− b − b ), ... , X i,i : ( b i −· · ·− b − b −· · ·− b i )(Fig.6),the molecules in the group B i have the edge monomers b i ( i = 1 , · · · , N ), and i is the biggest monomerindex number among the monomers included in the B i group molecule.In contrast to this, we classified network complexes into group LU j or group LD j according to theedge monomer index number, which is not the biggest index number in the molecule ( j = 1 , ..., N ) .Group LU j consists of molecules: X ,j : ( b j ), X ,j +1 : ( b j + 1 − b j ), ..., X N − j +1 ,N : ( b N − · · · − b j )(Fig.7). The molecules in the LU j group have monomer b j on the edge and no b − b polymer.Group LD j consists of molecules: X j,j : ( b j · · · b − b · · · b j ), X j +1 ,j +1 : ( b j +1 b j · · · b − b · · · b j ), ..., X N + j,N : ( b N · · · b j · · · b − b · · · b j ). The molecules in the LD j group have monomer b j on the edge and b − b polymer. According to the law of mass conservation, we determined N equations for relationships between the LU m +1 and ( B m + LD m ) reaction groups. First, we found a total concentration of monomer b N in the B ・・・ B m ・・・ B N X X ・・・ X ・・・ X X X ・・・ ・・・ ・・・ X
2, N ・・・ X ・・・ ・・・ ・・・ ・・・ ・・・ X ・・・ Xm,m ・・・ ・・・
N-1 ・・・ ・・・ ・・・ X N-1,N N ・・・ ・・・ ・・・ X N, N
N+1 ・・・ ・・・ X N+1, N ・・・
X2m,m ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ X X X Figure 4. Diagram of molecules in the N monomer network with the X l,m expression. The molecule X l,m has l monomersand the biggest monomer index number is m . When m ≥ l , the formation of molecular X l,m is ( b m − · · · − b m − l +1 ). When m < l , the formation of molecular X l,m is ( b m − · · · − b − b − · · · − b l − m ). LD m Group X(m+1, 2m+1) B m Group X(m, m)
X(m+i, m+i)
X(m+i, 2(m+i)) LU m+1 Group Figure 5. Reaction groups. The molecules in the ( b N − · · · − b − b − · · · − b N ) type network are classified into severalreaction groups. network, which is equal to initial concentration of [ b N (0)]. [ b N (0)] is derived as follows:[ b N (0)] = Σ Ni =1 X i,N ( t ) + X N,N ( t ) . (4)In the equation (4), X N,N ( t ) is doubled, because the complex, X(2N,N) has two b N monomers. [ b N (0)]is can also be written as follows: [ b N (0)] = Σ total B N + Σ total LD N . (5)Σ total B N is a total concentration of complexes in group B N , while Σ total LD N is the total concentrationof complexes in group LD N (Fig. 8).Following this, we derived the concentration of monomer [ b N − (0)]. In this network, the moleculesthat have b N − monomers, are included in the B N group and B N − group. Each molecule in B N − groupcontains b N − monomers. Each molecule in LD N − group contains two b N − monomers. No moleculesin LU N group have any b N − monomers. Therefore, we obtained the concentration of monomer [ b N − (0)] B ・・・ B m ・・・ B N-2 B N-1 B N X(1,1) X(1,2) ・・・
X(1,m) ・・・ ) X(1,N-2) X(1,N-1)
X(1,N) X(2,1) X(2,2) ・・・
X(2,m) ・・・ ・・・
X(2,N-1)
X(2, N) X(3,2) ・・・ ・・・ ・・・ ・・・ ・・・ ・・・
X(4,2) ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・
X(N-1,N)
X(2m, m) ・・・ ・・・ ・・・
X(N, N) ・・・ ・・・ ・・・
X(N+1, N) ・・・ ・・・ ・・・ ・・・ ・・・ ・・・
X(N+m, N) ・・・
X(2N-3,N-1) ・・・
X(2(N-1),N-1)
X(2N-2,N)
X(2N-1,N)
X(2N,N) B m Group B B N Group B N-1 Group B Group
Figure 6. B groups. The molecules in the ( b N − · · · − b − b − · · · − b N ) type network are classified to N B groups accordingto the biggest monomer index number of the complex. For example, B m group has 2 m molecules, which have ( b m −· · ·− b i )forms (monomer index i is smaller or equal to m ). by the addition of Σ total B N − , Σ total LD N − , and − Σ total LU N to the concentration of monomer [ b N (0)],as follows: [ b N − (0)] − [ b N (0)] = Σ total B N − + Σ total LD N − − Σ total LU N . (6)If the initial concentration [ b m +1 (0)] of b m +1 molecule is already known, we can calculate the initialconcentration [ b m (0)] of b m molecule as follows:[ b m (0)] − [ b m +1 (0)] = Σ total B m + Σ total LD m − Σ total LU m +1 . (7)Therefore, we were able to obtain N mass conservation relationships between the reaction groups. Wedetermine the sum of the equation (7), from m to N . Left-side terms of this equation cancel each other,leaving only [ b m (0)] − [ b N (0)]:[ b m (0)] − [ b N (0)] = Σ Nj = m ([ b j (0)] − [ b j +1 (0)]) , = Σ N − j = m (Σ total B j + Σ total LD j − Σ total LU j +1 ) . Followed with: [ b m (0)] = Σ Nj = m (Σ total B j + Σ total LD j ) − Σ N − j = m Σ total LU j +1 . (8) As previously described, mass preservation relationships between the reaction groups, B , LU , and LD groups are present. Here, we introduce the law of mass action for the reaction groups and their groupconcentrations at previous work. We show how to derive the mass action laws and their solutions in theappendix A.We define the parameter ξ m ( t ) as the total concentration of the molecules belonging to the LU m group. The parameter η m ( t ) is defined as the total concentration of molecules in B m and LD m groups. (cid:40) ξ m ( t ) = Σ total LU m ( t ) ,η m ( t ) = Σ total ( B m ( t ) + LD m ( t )) . (9) B ・・・ B m ・・・ B N-2 B N-1 B N X(1,1) X(1,2) ・・・
X(1,m) ・・・ ) X(1,N-2) X(1,N-1)
X(1,N)
X(2,1) X(2,2) ・・・
X(2,m) ・・・ ・・・
X(2,N-1)
X(2, N) X(3,2) ・・・ ・・・ ・・・ ・・・ ・・・ ・・・
X(4,2) ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・
X(N-1,N)
X(2m, m) ・・・ ・・・ ・・・
X(N, N) ・・・ ・・・ ・・・
X(N+1, N) ・・・ ・・・ ・・・ ・・・ ・・・ ・・・
X(N+m, N)
X(2N-3,N-1) ・・・
X(2(N-1),N-1)
X(2N-2,N)
X(2N-1,N)
X(2N,N) LD m Group LU m Group LU LD LU N Group LU N-1 Group LD N Group LD N-1 Group
Complexes in Lm group have monomer bm on right side edge. Figure 7. The molecules in the ( b N − · · · − b − b − · · · − b N ) type network are classified to N LU groups and N LDgroups, according to the edge monomer index number in the molecule, which is not the biggest monomer index number inthe molecule. LU m group has m molecules, which have b m monomer on the edge and no b − b dimer. Their molecularform is ( b i − · · · − b m ). The monomer index number i is bigger than or equal to m . LD m group has m molecules, whichhave b m monomer on the edge and b − b dimer form. The molecular form is ( b i − · · · − b − b − · · · − b m ) (monomerindex i is smaller or equal to m ). LD m Group B m Group LU m+1 Group Figure 8. N equations define the relationship between the LU m +1 and ( B m + LD m ) reaction groups, according to thelaw of mass conservation. The molecules in the blue region contain b m +1 monomers. The molecules in the purple regioncontain b m monomers. Afterward, equations (7)-(8) are rewritten using ξ m +1 ( t ) and η m ( t ):[ b m (0)] − [ b m +1 (0)] = η m ( t ) − ξ m +1 ( t ) . (10)The LU reaction group concentrations are summarized as follows: ξ N ( t ) = (cid:40) ξ ∗ N + − C N ξ ∗ n − e − βnt − C n e − βnt ( l N > , m = N ) , [ b N (0)] k N [ b N (0)] t +1 , ( l N = 0 and [ b N − (0)] = [ b N (0)] , m = N ) . (11) ξ m ( t ) = (cid:40) ξ ∗ m + − C m ξ ∗ m − e − βmt − C m e − βmt ( l m > , m : 1 ↓ N − , [ b m (0)] k m [ b m (0)] t +1 ( l m = 0 and [ b m − (0)] = [ b m (0)] , m : 1 ↓ N − . (12) ( t ) = ξ ∗ − C ξ ∗ − e − β t − C e − β t ( l > , m = 1) , [ b (0)] k [ b (0)] t +1 ( l = 0 , m = 1) . (13)LU group concentrations, ξ N , ..., ξ m , ..., ξ , are independent of each other and have strict solutions, whichconverge to equilibrium solutions when time t tends to + ∞ .The total concentration of B m and LD m groups η m is derived from the equation (7), as follows: (cid:40) η N ( t ) = [ b N (0)] , ( m = N ) ,η m ( t ) = [ b m (0)] − [ b m +1 (0)] + ξ m +1 ( t ) , ( m = 1 , · · · , N − . (14)Here, we showed that the total concentration of LU m group, or total concentrations of LD m and B m groups, can be solved. They converge to the equilibrium solutions when time t tends to + ∞ .We regarded the reaction group concentrations as a priori upper boundaries of the concentration ofeach molecular in the group. Afterward, we demonstrated that each molecule in a group converges to astable solution. All evolution equations of the complexes in the pathway network were derived. We show that ODEsof all concentrations of the molecules X l,m ( t ) in the network have the following shape: dX l,m ( t ) dt = − A l,m ( t ) X l,m ( t ) + f l,m ( t ) . (15)First, we solved the homogeneous differential equation, dX l,m ( t ) dt = − A l,m ( t ) X l,m ( t ), and obtained thesolution X l,m ( t ) = X l,m (0) e − (cid:82) t A l,m ( s ) ds . These equations were solved by varying the parameters. X l,m ( t ) = e − (cid:82) t A l,m ( s ) ds + (cid:90) t e − (cid:82) s A l,m ( u ) du f l,m ( s ) ds. (16)Afterward, all complex concentrations were shown to be integrable.We set equilibrium values for molecules X l,m , A m ( t ), and f m ( t ) as X ∗ l,m , A ∗ l,m , and f ∗ l,m . We definedthe parameter Y ( t ) = X l,m ( t ) − X ∗ l,m : dY l,m ( t ) dt = − A l,m ( t ) Y l,m ( t ) + (cid:32) f l,m ( t ) − f ∗ l,m A ∗ l,m A l,m ( t ) (cid:33) . (17)We defined g l,m ( t ) = f l,m ( t ) − f ∗ l,m A ∗ l,m A l,m ( t ) and a l,m ( t ) = (cid:82) tT A l,m ( s ) ds . Y l,m ( t ) = e − a ( t )+ a ( T ) Y l,m ( T ) + (cid:90) tT e a ( s ) − a ( t ) g l,m ( s ) ds ≤ e a ( T ) | Y l,m ( T ) | e − A ∗ l,m ( t − T ) + (cid:90) tT e − A ∗ l,m ( s − T ) e − at e − δs ds, (18)where a ( t ) ≥ A ∗ l,m ( t − T ).Because A l,m ( t ) = k i ξ i ( t ) + k j ξ j ( t ) + l i + l j and integral of ξ i ( t ) is written as follows: (cid:90) tT ξ m ( s ) ds = (cid:20) ξ ∗ m + s + 1 B m log (1 − C m ) e − B m s (cid:21) tT = ξ ∗ m + ( t − T ) + 1 B m log (cid:18) − C m e − B m t − C m e − B m T (cid:19) , (19)where B m = k m ( ξ ∗ m + − ξ ∗ m − ), C m = [ b m (0)] − ξ ∗ m + [ b m (0)] − ξ ∗ m − . When time t tends to ∞ , Y l,m ( t ) tends to 0:lim t →∞ X l,m ( t ) = f ∗ l,m A ∗ l,m . (20) .1 The concentration and stability of complexes in LU groups We calculated the concentrations of the molecules by starting with the complex X , ⊂ LD , B followedby solving ODEs for complexes X , ⊂ B , X , ⊂ B , · · · ., X N − m +1 ,N ⊂ B N in LD group, step by step.Afterward, we solved the ODEs of complexes in LD , · · · , LD N group.Additionally, we began with the molecules in the LU N group, and solved the concentrations ofmolecules from LU N − , · · · , LU , LD , · · · , LD N groups, step by step. LU N group First, we started with LU N group and calculated the complex concentrations of LU N − , · · · , LU , LU . LU N group has only one complex, X ,N , and its concentration X ,N ( t ) is equal to LU N group solution, ξ N ( t ). X ,N ( t ) = ξ N ( t ) . (21)Its equilibrium is obtained. X ∗ ,N = ξ ∗ N . (22) LU N − ( N ≥ ) group We considered the concentrations of molecules in LU N − ( N ≥
2) group, X i,N + i − ( t ) ⊂ B i (i = 1, 2 ). dX i,N + i − ( t ) dt = − A i,N + i − ( t ) X i,N + i − ( t ) + f i,N + i − ( t ) ( i = 1 ,
2) (23)where A i,N + i − ( t ) = k N ξ N ( t ) + k N − η N − ( t ) + l N + l N − ( i = 1 , .f i,N + i − ( t ) = l N ξ N − ( t ) + l N − η N − ( t ) ( i = 1) ,k N ξ N ( t ) ξ N − ( t ) + l N − ( η N ( t ) − ξ N ( t )) ( i = 2) . (24)Its equilibrium solution is lim t →∞ X ∗ i,m + i − = f ∗ i,m + i − A ∗ i,m + i − ( i = 1 , , (25)where A ∗ i,N + i − = k N ξ ∗ N + k N − η ∗ N − + l N + l N − ( i = 1 , ,f ∗ i,N + i − = l N ξ ∗ N − + l N − η ∗ N − ( i = 1) ,k N ξ ∗ N ξ ∗ N − + l N − ( η ∗ N − ξ ∗ N ) ( i = 2) . (26) LU m group ( m = 2 , · · · , N − Afterward, we considered the concentrations of molecules in LU m group, specifically, of the LU m groupcomplex, X i,m + i − ( t ) ⊂ B i ( i = 1 , · · · , N − m + 1). Starting from the mass conservation laws and massaction laws, we obtained the reaction ODEs for X i,m + i − molecule as follows: dX i,m + i − ( t ) dt = − A i,m + i − ( t ) X i,m + i − ( t ) + f i,m + i − ( t )( i = 1 , · · · , N − m + 1) (27)here A i,m + i − ( t ) = k m +1 ξ m +1 ( t ) + k m η m − ( t ) + l m +1 + l m ( i = 1) ,k m + i ξ m + i ( t ) + k m η m − ( t ) + (cid:80) ij =0 l m + j ( i = 2 , · · · , N − m ) ,k m η m − ( t ) + (cid:80) N − mj =0 l m + j ( i = N − m + 1) ,f i,m + i − ( t ) = l m +1 ξ m ( t ) + l m η m ( t ) ( i = 1) ,l m + i ( ξ m ( t ) − (cid:80) i − j =1 X j,m + j − ( t )) + l m ( η m + i ( t ) − (cid:80) i − j =1 X j,m + i − ( t ))+ (cid:80) i − j =1 k m + j X i − j,m + i − ( t ) X j,m + j − ( t ) ( i = 2 , · · · , N − m ) , (cid:80) N − mj =1 k m + j X j,m + j − ( t ) X N − m +1 − j,N ( t ) + l m ( η N ( t ) − (cid:80) N − mj =1 X j,N ( t ))( i = N − m + 1) . Its equilibrium solution islim t →∞ X i,m + i − ( t ) = f ∗ i,m + i − A ∗ i,m + i − ( i = 1 , · · · , N − m + 1) , (28)where A ∗ i,m + i − = k m +1 ξ ∗ m +1 + k m η ∗ m + l m +1 + l m ( i = 1) ,k m + i ξ ∗ m + i + k m η ∗ m − + (cid:80) m + ij = m l j ( i = 2 , · · · , N − m ) ,k m η ∗ m − + (cid:80) N − mj =0 l m + j ( i = N − m + 1) ,f ∗ i,m + i − = l m +1 ξ ∗ m + l m η ∗ m ( i = 1) ,l m ( η ∗ m + i − (cid:80) i − j =1 X ∗ j,m + i − ) + l m + i ( ξ ∗ m − (cid:80) i − j =1 X ∗ j,m + j − )+ (cid:80) i − j =1 k m + j X ∗ i − j,m + i − X ∗ j,m + j − ( i = 2 , · · · , N − m − , (cid:80) N − mj =1 k m + j X ∗ j,m + j − X ∗ N − m − j +1 ,N + l m ( η ∗ N − (cid:80) N − mj =1 X ∗ j,N )( i = N − m + 1) . LU group complexes Reaction ODE of LU complexes, X i,i ⊂ B i ( i = 1 , · · · , N ), is written as follows: dX , ( t ) dt = − A , ( t ) X , ( t ) + f , ( t ) , (29)where (cid:40) A , ( t ) = 2 k ξ ( t ) + 2 k ξ ( t ) + (cid:80) j =1 l j ,f , ( t ) = l η ( t ) + l ξ ( t ) . Therefore, the equilibrium solution of X i,i ( t ) is X ∗ , = lim t →∞ X , ( t ) = f ∗ , A ∗ , , (30)here (cid:40) A ∗ , = 2 k ξ ∗ + 2 k ξ ∗ + (cid:80) j =1 l j ,f ∗ , = l η ∗ + l ξ ∗ . (31)When i ¿ 1, dX i,i ( t ) dt = − A i,i ( t ) X i,i ( t ) + f i,i ( t ) , ( i = 2 , · · · , N − , (32)where A i,i ( t ) = 2 k ξ ( t ) + k i +1 ξ i +1 ( t ) + (cid:80) i +1 j =1 l j ,f i,i ( t ) = (cid:80) i − j =1 k j +1 X j,j ( t ) X i − j,i ( t ) + l i +1 { ξ ( t ) − (cid:80) i − j =1 X j,j ( t ) } + l { [ b i (0)] − [ b i +1 (0)] + ξ i +1 ( t ) − (cid:80) i − j =1 X j,i ( t ) } . Therefore, the equilibrium solution of X i,i ( t ) is X ∗ i,i = lim t →∞ X i,i ( t ) = f ∗ i,i A ∗ i,i , (33)where A ∗ i,i = 2 k ξ ∗ + k i +1 ξ ∗ i +1 + (cid:80) i +1 j =1 l j ,f ∗ i,i = (cid:80) i − j =1 k j +1 X ∗ j,j X ∗ i − j,i + l i +1 { ξ ∗ − (cid:80) i − j =1 X ∗ j,j } + l { [ b i (0)] − [ b i +1 (0)] + ξ ∗ i +1 − (cid:80) i − j =1 X ∗ j,i } . (34)Finally, for the molecules belonging to the LU , X N,N ⊂ B N , reaction ODE of this molecule is written: dX N,N ( t ) dt = − A N,N ( t ) X N,N ( t ) + f N,N ( t ) , (35)where (cid:40) A N,N ( t ) = 2 k ξ ( t ) + (cid:80) nj =1 l j ,f N,N ( t ) = (cid:80) N − j =1 k j +1 X j,j ( t ) X N − j,N ( t ) + l ([ f N (0)] − (cid:80) N − j =1 X j,N ( t )) . Therefore, the equilibrium solution of X ∗ N,N = lim t →∞ X N,N ( t ) is X ∗ N,N = f ∗ N,N A ∗ N,N , (36)where A ∗ N,N = 2 k ξ ∗ + (cid:80) Nj =1 l j ,f ∗ N,N = (cid:80) N − j =1 k j +1 X j,j ( t ) X ∗ N − j,N + l ([ b N (0)] − (cid:80) N − j =1 X ∗ j,N ) . We demonstrated that the complex concentrations in LU groups are integrable and converge to theequilibrium. Following this, we solved the concentrations of LD group complexes. Here, we initiallydetermined the concentrations of X , ⊂ LD , B complex, followed by ODEs for complexes X , ⊂ B , X , ⊂ B , · · · , X N − m +1 ,N ⊂ B N in LD group. Finally, we solved the ODEs of complexes belongingto LD group, · · · , LD N group, step by step. .7 Concentration of complexes in LD group LD group contains N complexes, X , , X , , · · · , X N +1 ,N . This is a group with special characteristics,as the complexes belonging to this group have b − b connection on the edge. Any complex in LD canreact with every complex in LU group.First, we considered the concentration of LD group complex, X i +1 ,i ( t ) ⊂ B i ( i = 1 , · · · , N ). Fromthe law of mass conservation and the law of mass action, we obtained reaction ODEs of molecule X i +1 ,i as follows: dX i +1 ,i ( t ) dt = − A i +1 ,i ( t ) X i +1 ,i ( t ) + f i +1 ,i ( t ) ( i = 1 , · · · , N ) , (37)where A i +1 ,i ( t ) = k ξ ( t ) + 2 l + l ( i = 1) ,k ξ ( t ) + k i +1 ξ i +1 ( t ) + l + 2 l + (cid:80) i +1 j =3 l j ( i = 2 , · · · , N − ,k ξ ( t ) + l + 2 l + (cid:80) Nj =3 l j ( i = N ) ,f i +1 ,i ( t ) = k X , ( t ) + l ( η ( t ) − X , ( t )) ( i = 1) , k X , ( t ) X i,i ( t ) + k X , ( t ) X i − ,i ( t ) + (cid:80) i − j =1 k j +1 X j +1 ,j ( t ) X i − j,i ( t )+ l { η i ( t ) − (cid:80) ij =1 X j,i ( t ) } + l i +1 ( η ( t ) − (cid:80) j =1 X j, ( t ) − (cid:80) i − j =1 X j +1 ,j ( t ))( i = 2 , · · · , N − , k X , ( t ) X N,N ( t ) + 2 k X , X i − ,i + (cid:80) N − j =2 k j +1 X j +1 ,j ( t ) X N − j,N ( t )+ l [ b N (0)] − (cid:80) Ni =1 X i,N ( t ) ( i = N ) . Deformation of equation (37) is: ddt { X i +1 ,i ( t ) − b ∗ i +1 ,i A ∗ i +1 ,i } = − A i +1 ,i ( t ) { X i +1 ,i ( t ) − b ∗ i +1 ,i A ∗ i +1 ,i } + f i +1 ,i ( t ) − b ∗ i +1 ,i A ∗ i +1 ,i A i +1 ,i ( t )( i = 1 , · · · , N ) . (38)Let Y i +1 ,i ( t ) be a difference between X i +1 ,i ( t ) and its equilibrium solution, X ∗ i +1 ,i = f ∗ i +1 ,i A ∗ i +1 ,i . Y i +1 ,i ( t ) isequal to X i +1 ,i ( t ) − X ∗ i +1 ,i .Then, ddt Y i +1 ,i ( t ) = − A i +1 ,i ( t ) Y i +1 ,i ( t ) + f i +1 ,i ( t ) − f ∗ i +1 ,i A ∗ i +1 ,i A i +1 ,i ( t ) . (39)General and equilibrium solutions were obtained as follows: Y i +1 ,i ( t ) = e − a ( t )+ a ( T ) Y i +1 ,i ( T ) ( f i +1 ,i ( t ) A i +1 ,i ( t ) = f ∗ i +1 ,i A ∗ i +1 ,i ) ,e − a ( t )+ a ( T ) Y i +1 ,i ( T ) + e − a ( t ) (cid:82) ts = T e a ( s ) g ( s ) ds ( a i +1 ,i ( t ) dt = A i +1 ,i ( t ) , g i +1 ,i ( t ) = f i +1 ,i ( t ) − f ∗ i +1 ,i A ∗ i +1 ,i A i +1 ,i ( t )) . (40)Therefore equilibrium solutions of X ∗ i +1 ,i = lim t →∞ X i +1 ,i ( t ) are: X ∗ i +1 ,i = f ∗ i +1 ,i A ∗ i +1 ,i , (41)here A ∗ i +1 ,i = k ξ ∗ + 2 l + l ( i = 1) ,k ξ ∗ + k i +1 ξ ∗ i +1 + l + 2 l + (cid:80) i +1 j =3 l j ( i = 2 , · · · , N − ,k ξ ∗ + l + 2 l + (cid:80) nj =3 l j ( i = N ) ,f ∗ i +1 ,i = k X ∗ , + l ( η ∗ − X ∗ , ) ( i = 1) , k X ∗ , X ∗ i,i + k X ∗ , X ∗ i − j,i + (cid:80) i − j =1 k j +1 X ∗ j +1 ,j X ∗ i − j,i + l ( η ∗ i − (cid:80) ij =1 X ∗ j,i )+ l i +1 ( η ∗ − (cid:80) j =1 X ∗ j, − (cid:80) i − j =1 X ∗ j +1 ,j ) ( i = 2 , · · · , N − , k X ∗ , X ∗ N,N + (cid:80) N − j =1 k j +1 X ∗ j +1 ,j X ∗ N − j,N + l { [ b N (0)] − (cid:80) Ni =1 X ∗ i,N } ( i = N ) . LD m group LD m ( m = 2 , · · · , N −
1) group has N − m + 1 complexes, X m,m , · · · ,X m + i,m + i , · · · , X N + m +1 ,N . We solved the concentrations of complex from m = 2 to N −
1. A complexin LD m group can react with every complex in LU m +1 group.From the mass conservation law and mass action law, we obtain reaction ODEs for the concentrations of LD m group complex, X m + i,m + i ( t ) ⊂ B i . dX m + i,m + i ( t ) dt = − A m + i,m + i ( t ) X m + i,m + i ( t ) + f m + i,m + i ( t ) ( i = 0 , · · · , N ) . (42)First, we considered the concentration of the complex, X m,m ( t ) ⊂ B m , LD m (Fig.9). This complex issymmetric, with b m monomers on both edges, b m . . . b b . . . b m . The reactions between the complexesin B N group lead to the generation of X m,m complex. The dissolution of LD m group complexes (e.g., X m + i,m + i ( i = 1 , · · · , N − m )), leads to the generation of X m,m complex and LU m +1 group complexes(e.g., X m +1+ i,i ( i = 1 , · · · , N − m − LD m Group X(m+1, 2m+1) B m Group X(m, m)
X(m+i, m+i)
X(m+i, 2(m+i))
Figure 9. Groups related to the production of complex X m,m ⊂ LD m in the network diagram. Molecules in the greenregion of B m group react each other and produce complex X m,m . Molecules in the pink region of LD m group are dissolvedand complex X m,m is produced. he coefficients of the equation (42), A m,m and b m,m are defined as follows: A m,m ( t ) = 2 k m +1 ξ m +1 ( t ) + 2 (cid:80) m +1 j =1 l j ( i = 0) ,f m,m ( t ) = 2 k X m,m ( t ) + (cid:80) m − j =1 k j +1 X m − j,m ( t ) X m + j,m ( t )+ l m +1 ( η m ( t ) − (cid:80) m − j =1 X j,m ( t )) ( i = 0) . (43)Additionally, we considered the concentration of the complex, X m + i,m + i ⊂ B m + i ( i = 2 , · · · , N − m −
1) (Fig.4 . LD m Group X2m+1, m+1 B m Group Xm, m X m+i, m+i X2(m+i), m+i B m+i Group LD m+i Group X2m+i, m+i Figure 10. Groups related to the production of complex X m + i,m + i ⊂ LD m in the network diagram. Molecules in the greenregion of B m and LD m groups and molecules in the green region of B m + i group react each other and produce complex X m,m . Molecules in the pink region of LD m , B m + i and LD m + i groups are dissolved and produce complex X m + i,m + i . A m + i,m + i ( t ) = k m +1 ξ m +1 ( t ) + k m + i +1 ξ m + i +1 ( t ) + l + 2 (cid:80) m +1 j =2 l j + (cid:80) m + i +1 j = m +2 l j ,f m + i,m + i ( t ) = 2 k X m + i,m + i ( t ) X m,m ( t ) + l m + i ( η m + i ( t ) − (cid:80) m + i − j =1 X j,m + i ( t ))+ (cid:80) m +1 j =2 k j X m + i − j +1 ,m + i ( t ) X m + j − ,m ( t )+ (cid:80) m + ij = m +1 k j X m + i − j +1 ,m + i ( t ) X m + j − ,j − ( t )+ (cid:80) mj =2 k j X m + i + j − ,m + i ( t ) X m − j +1 ,m ( t )+ l m +2 ( η m ( t ) − (cid:80) mj =1 X j,m ( t ) − (cid:80) i − j =0 X m + j,m + j ( t )) . (44)IIf we consider the concentration of the complex X N + m,N ( t ) ⊂ LDm , we observe that the X N + m,N complex is produced as a consequence of a reaction between a complex belonging to the B m group anda complex belonging to the B N group, or between a complex in the LD m group and another one in the B N group, and by degradation of the bigger complexes in B N group(Fig.11)Therefore, the coefficients of equation (42) were derived as follows: A N + m,N ( t ) = k m +1 ξ m +1 ( t ) + l + 2 (cid:80) m +1 j =2 l j + (cid:80) Nj = m +2 l j ( m < N + 1) ,k m +1 ξ m +1 ( t ) + l + 2 (cid:80) m +1 j =2 l j ( m = N + 1) ,f N + m,N ( t ) = (cid:40) (cid:80) m − j =0 k j +1 X m − j,m ( t ) X N + j,N ( t ) + sum mj =0 k j +1 X m + j,m ( t ) X N − j,N ( t )+ (cid:80) N − j = m k j +1 X m + j,j ( t ) X N − j,N ( t ) + l m +1 ([ b N (0)] − (cid:80) N + m − j =1 X j,N ) . D m Group X(m+1, 2m+1) B m Group X(m, m)
X(m+i, m+i) B N Group
Figure 11. Groups related to the generation of complex X N + m,N ⊂ LD m in the network diagram. Molecules in the greenregion of B m and LD m groups and molecules in the green region of B N group react each other and produce complex X m,m . Molecules in the pink region of LD N and B N groups are dissolved and produce complex X N + m,N . We have already obtained the values and equilibrium solutions of the parameters ξ m +1 ( t ) , ξ m + i +1 ( t ) , X j,m ( t )( j =1 , · · · , m ), X m + i,m + i ( t )( i = 1 , · · · , N − m ). It is possible to solve the equation (42), step by step, from i = 0 to N, as previously done, and obtain the equilibrium solutions. X ∗ m + i,m + i = b ∗ m + i,m + i A ∗ m + i,m + i ( i = 1 , · · · , N − m ) , (45)where A ∗ m + i,m + i = k m +1 ξ ∗ m +1 + 2 (cid:80) m +1 j =1 l j ( i = 0) ,k m +1 ξ ∗ m +1 + k m + i +1 ξ ∗ m + i +1 + l + 2 (cid:80) mj =2 l j + (cid:80) m + ij = m +1 l j ( i = 1 , · · · , N − m − ,k m +1 ξ ∗ m +1 + l + 2 (cid:80) m +1 j =2 l j + (cid:80) Nj = m +2 l j ( i = N − m, m < N − ,k m +1 ξ ∗ m +1 + l + 2 (cid:80) m +1 j =2 l j ( i = N − m, m + 1 = N ) (46) f ∗ m + i,m + i = k X ∗ m,m + (cid:80) m − j =1 k j +1 X ∗ m − j,m X ∗ m + j,m + l m +1 ( η ∗ m − (cid:80) m − j =1 X ∗ j,m )( i = 0) , k X ∗ m + i,m + i X ∗ m,m + (cid:80) m +1 j =2 k j X ∗ m + i − j +1 ,m + i X ∗ m + j − ,m + (cid:80) m + ij = m +1 k j X ∗ m + i − j +1 ,m + i X ∗ m + j − ,j − + (cid:80) mj =2 k j X ∗ m + i + j − ,m + i X ∗ m − j +1 ,m + l m +2 ( η ∗ m − (cid:80) mj =1 X ∗ j,m − (cid:80) i − j =0 X ∗ m + j,m + j )+ l m + i ( η ∗ m + i − (cid:80) m + i − j =1 X ∗ j,m + i )( i = 2 , · · · , N − m − , (cid:80) m − j =0 k j +1 X ∗ m − j,m X ∗ N + j,N + (cid:80) mj =0 k j +1 X ∗ m + j,m X ∗ N − j,N + (cid:80) N − j = m k j +1 X ∗ m + j,j X ∗ N − j,N + l m +1 ([ b N (0)] − (cid:80) N + m − j =1 X ∗ j,N )( i = N − m ) . .9 Concentration of complexes in LD N group LD N group consists of only one complex, X N,N . This complex can be produced only in the reactionsbetween complexes in B N group, but it cannot be obtained as a product of degradation, because itrepresents the biggest complex in the network.From mass conservation law and mass action law, we obtained the reaction ODE of the X N,N complexas follows: dX N,N ( t ) dt = − A N,N ( t ) X N,N ( t ) + f N,N ( t ) , (47)where (cid:40) A N,N ( t ) = l + 2 (cid:80) Nj =2 l j ,f N,N ( t ) = 2 k X N,N ( t ) + (cid:80) N − j =1 k j +1 X j,N ( t ) X N − j,N ( t ) . Coefficient A N,N ( t ) is constant and positive, as shown in the equation (48). We set this coefficient as A ∗ N,N = l + 2 (cid:80) Nj =2 l j , which allowed us to solve the equation (47) by varying the parameters, and X N,N (0) = 0. X N,N ( t ) = (cid:90) t e − ( t − s ) A ∗ N,N f N,N ( s ) ds. (48)When the parameter Y N,N ( t ) is set as Y N,N ( t ) = X N,N ( t ) − f ∗ N,N A ∗ N,N , from the equation (47), we canobtain: dY N,N ( t ) dt = − A ∗ N,N Y N,N ( t ) + f N,N ( t ) − f ∗ N,N . (49)We obtained the equilibrium solutions as previously: X ∗ N,N = f ∗ N,N A ∗ N,N , ( i = 1 , · · · , N − m ) , (50)where (cid:40) A ∗ N,N = l + 2 (cid:80) Nj =2 l j f ∗ N,N = 2 k X ∗ N,N + (cid:80) n − j =1 k j +1 X ∗ j,N X ∗ N − j,N . This led to a conclusion that all complexes in the ( b N −· · ·− b − b −· · ·− b N ) type network are integrableand have explicit equilibrium solutions. We applied the equilibrium values of complex concentrations, determined for the N monomer network,to the 3-monomer MMP2/TIMP2/MT1-MMP network. The equilibrium values of the group and complexconcentrations of this network were calculated using R ver. 3.2.2 wuth desolve package.The parameter X , ( t ) is a concentration of the complex (MMP2/TIMP2/MT1-MMP/MT1-MMP)at time t , and b corresponds to MT1-MMP, b to TIMP2, b to MMP2.From this simulation, we can see that the theoretical result of the target complex concentration X , ( ∞ ) agrees with the results obtained by the ODE simulations of the target complex concentration X , ( ∞ ) (Fig.12). The group concentration η ( ∞ ) was shown to be always bigger than X , ( ∞ ).The difference between monomers in the network, [ b (0)] − [ b (0)], represents a regulatory parameterfor the η group, ξ group, and the complexes that belong to these groups, such as X , . When [ b (0)] − [ b (0)] < η ( ∞ ) and X , ( ∞ ) are positive and ξ ( ∞ ) = 0. When [ b (0)] − [ b (0)] = 0, η ( ∞ ), ξ ( ∞ ),and X , ( ∞ ) are 0. When [ b (0)] − [ b (0)] > η ( ∞ ) = X , ( ∞ ) = 0, and ξ ( ∞ ) is positive. 𝝃 (∞) 𝜼 (∞) 𝑿 (∞) Mathematical [b1(0)]-[b2(0)] 𝑿 (∞) Simulation [b1(0)]-[b2(0)]=0 Figure 12. Simulation and theoretical results of the equilibrium concentrations of the molecule X , , η group, and ξ groupat [ b (0)] = 1000 nM , during the MMP2 activation process. X , corresponds to the (MMP2/TIMP2/MT1-MMP/MT1-MMP) complex. The black dotted line represents X , ( ∞ ), according to ODE simulations. Red line represents theconcentration of X , ( ∞ ) molecule, calculated using the equation(68). The results we obtained and the simulation resultshow a high level of agreement. Green line represents the equilibrium group concentration η ( ∞ ), where the theoreticalequations were applied. η group includes the target molecule X , . Blue line represents the equilibrium group concentration ξ ( ∞ ), where we applied our theoretical equations. ξ group has dual relationship with η group. The positive or negativedifference in the concentrations, [ b (0)] − [ b (0)], determines the positive and zero η ( ∞ ) and ξ ( ∞ ) groups. This demonstrates that this is crucial for the determination of η ( ∞ ), ξ ( ∞ ), and X , ( ∞ ) concen-trations.Figure 13 shows X , ( ∞ ) when [ b (0)] and [ b (0)] vary from 0 to 1000 nM. For [ b (0)] − [ b (0)] > X , ( ∞ ) is always 0. By generalizing MMP2 activation network to N monomer network with ( b N · · · b − b · · · b N ) com-plexes, we can observe the group reaction behavior and the layered structure of the network. This allowedus to take a two-step approach, in order to solve the pathway network system.As the first step, we classified the complexes to the reaction groups, such as B m , LU m , and LD m ,and these groups behave according to the mass conservation laws. Therefore, we aggregated complexODEs into the reaction group ODEs, and obtained the strict solutions for reaction group concentrations.In the following step, we solved all complex concentrations in the each reaction group, while the groupconcentrations were obtained at the first step. The concentrations of all complexes in the N monomerpathway network are derived explicitly and converge to equilibria.We think our approach is useful for biologists dealing with the pathway network with homodimersymmetric formation complexes. We show this type pathway network is integrable and stable. Eachconcentration of the complexes in the network is solved explicitly.Especially, we obtain equilibrium con-centration of each complex easily and quickly with our approach.We also identify the difference between the initial concentrations of the reacting monomers, [ b m (0)] − [ b m +1 (0)], as regulator parameter of the target complex behavior. For example, in the case of MMP2activation pathway network, the difference of the initial concentration between monomer b :MT1-MMPand monomer b :TIMP2, [ b (0)] − [ b (0)] is regulator of the target complex X (4 , b2(0)] [b1(0)] X4,3(∞)
Figure 13. The result of the equilibrium concentration simulation for the molecule X , when [ b (0)] and [ b (0)] vary from0 nM to 1000 nM during the MMP2 activation process. When not all the reaction rate constants and initial concentrations but the relevant constants andinitial concentrations are known, we can calculate the group concentrations the target complex related,as the upper value of the target complex. It is important to obtain the upper value of the target complexconcentration in the difficult case to get all informations of the pathway network. It is difficult to simulatethe target complex concentration directly with the ODEs.Our mathematical approach may help us understand the mechanism of this type pathway network byknowing the background mathematical laws which govern the network. We think that pathway networksand chemical networks with homodimer symmetric complex formation ( b N · · · b b · · · b N ) have similardominant mathematical laws. Autophagy complex and its effects on the production of Atg-PE8 are verysimilar to MMP2 activation and (MT1-MMP/TIMP2/MMP2) complex formation. In this research, wetake account of only chemical reactions of monomers and complexes for simplicity. Competing interests
The authors declare that they have no competing interests.
Author’s contributions
The original mathematical models were drawn from earlier work of TS. Modifications were made byKI. Modeling and computational works were carried out by the first author with suggestions from thesecond. All other works were executed jointly.
Acknowledgements
This work is founded by JSPS Core-to-Core Program “Establishing International Research Networkof Mathematical Oncology”.
Group ConcentrationA.1 LU m group As we saw previously, the reaction groups follow mass preservation laws. LU m +1 group has dual rela-tionships with B m and LD m groups, as seen in equation (7) ( m = 1 , · · · , N − LU m +1 group as: d (cid:80) total LU m +1 ( t ) dt = N − m (cid:88) i =1 dX i,m + i ( t ) dt = − k m +1 (cid:88) total LU m +1 ( t ) { (cid:88) total B m ( t ) + (cid:88) total LD m ( t ) } − l m +1 (cid:88) total LU m +1 ( t )+ l m +1 { N − m (cid:88) j =1 ( (cid:88) total B m + j ( t ) + (cid:88) total LD m + j ( t )) − N − m (cid:88) j =2 (cid:88) total LU m + j ( t ) } . We define the parameter ξ m ( t ) as the total concentration of the molecules belonging to the LU m group.The parameter η m ( t ) is defined as the total concentration of molecules in B m and LD m groups. (cid:40) ξ m ( t ) = Σ total LU m ( t ) ,η m ( t ) = Σ total ( B m ( t ) + LD m ( t )) . (51)Afterward, equations (7)-(8) are rewritten using ξ m +1 ( t ) and η m ( t ):[ b m (0)] − [ b m +1 (0)] = η m ( t ) − ξ m +1 ( t ) . (52)[ b m +1 (0)] = Σ Nj = m +1 η j ( t ) − Σ N − j = m +1 ξ j +1 ( t ) . (53)Then, dξ m +1 ( t ) dt = − k m +1 ξ m +1 ( t ) − { ([ b m (0)] − [ b m +1 (0)]) k m +1 + l m +1 } ξ m +1 ( t )+ l m +1 [ b m +1 (0)] . (54)We obtain LU group ODEs in this way ( m = 1 , · · · , N ). N ( N + 1) molecule ODEs were aggregatedto N reaction group ODEs, which can be solved explicitly and have asymptotically stable solutions.The discriminant of dξ m +1 ( t ) dt is: D m +1 = { k m +1 ([ b m (0)] − [ b m +1 (0)]) + l m +1 } + 4 k m +1 l m +1 [ b m +1 (0)] . (55)If D m +1 is positive, there are two equilibrium solutions of dξ m +1 ( t ) dt = 0. ξ ∗ m +1 ± = −{ k m +1 ([ b m (0)] − [ b m +1 (0)]) + l m +1 } ± (cid:112) D m +1 k m +1 . (56)The explicit solution of equation(54) is: ξ m +1 ( t ) = ξ ∗ m +1+ − C m +1 ξ ∗ m +1 − e − β m +1 t − C m +1 e − β m +1 t , (57)where C m +1 = [ b m +1 (0)]+ ξ ∗ m +1 − [ b m +1 (0)]+ ξ ∗ m +1+ , β m +1 = k m +1 ( ξ ∗ m +1+ − ξ ∗ m +1 − ).hen time t tends to + ∞ , e − β m +1 t → +0.lim t →∞ ξ m +1 ( t ) = ξ ∗ m +1+ . (58)If l m +1 = 0 and [ b m (0)] = [ b m +1 (0)], the discriminant D m +1 = 0. In this case, we obtain a steadysolution of equation(54) as ξ ∗ m +1 = 0. The equation(54) is rewritten as follows: dξ m +1 ( t ) dt = − k m +1 ξ m +1 ( t ) . (59)The explicit solution of equation(54) in this case is: ξ m +1 ( t ) = [ b m +1 (0)] k m +1 [ b m +1 (0)] t + 1 . (60) ξ m +1 ( t ) decreases with O ( t − ) more slowly than in the case where D m +1 > O ( e − β m +1 t ). Whentime t tends to + ∞ , the numerator is constant and denomination tends to + ∞ :lim t →∞ ξ m +1 ( t ) = ξ ∗ m +1 = 0 . After obtaining LU m +1 group concentration ξ m +1 ( t ) and its equilibrium value ξ ∗ m +1 , we obtained the B m group and LD m group concentration η m ( t ) and its equilibrium value η ∗ m , using the mass conservationlaw equation (52). η m ( t ) = [ b m (0)] − [ b m +1 (0)] + ξ m +1 ( t ) . (61) η ∗ m = [ b m (0)] − [ b m +1 (0)] + ξ ∗ m +1 . (62) LU N and LU reaction group ODEs are somewhat different from LU m reaction groups ( m = 2 · · · , N − A.2 LU N Group LU N group consists of only molecule X (1 , N ), and the ODE of the molecule X (1 , N ) is: dX ,N ( t ) dt = − k N X ,N ( t )( (cid:88) total B N − + (cid:88) total LD N − )+ l N ( (cid:88) total B N + (cid:88) total LD N − X ,N ( t )) . (63)The LU N ( t ) reaction group concentration is equal to the molecule concentration X ,N ( t ). Therefore,the ODE of the LU N reaction group is: dξ N ( t ) dt = − k N ξ N ( t )([ b N − (0)] − [ b N (0)] + ξ N ( t )) + l N ([ b N (0)] − ξ N ( t ))= − k N ξ N ( t ) − { k N ([ b N − (0)] − [ b N (0)]) + l N } ξ N ( t ) + l N [ b N (0)] . (64)The solutions for the equation (64) are obtained in the similar way as LU m reaction group concentration.The discriminant of dξ N ( t ) dt is: D N = { k N ([ b N − (0)] − [ b N (0)]) + l N } + 4 k N l m +1 [ b N (0)] . (65)If D N is positive, there are two equilibrium solutions of dξ N ( t ) dt = 0. ξ ∗ N ± = −{ k N ([ b N − (0)] − [ b N (0)]) + l N } ± (cid:112) D N k N . (66)he explicit solution of equation (64) is: ξ N ( t ) = ξ ∗ N + − C N ξ ∗ N − e − β N t − C N e − β N t . (67)where C N = [ b N (0)]+ ξ ∗ N − [ b N (0)]+ ξ ∗ N + , β N = k N ( ξ ∗ N + − ξ ∗ N − ).When time t tends to + ∞ , e − β N t → +0: lim t →∞ ξ N ( t ) = ξ ∗ N + . (68)In the case when l N = 0 and [ b N − (0)] = [ b N (0)], i.e., the discriminant D N = 0, the equation (64) isrewritten as follows: dξ N ( t ) dt = − k N ξ N ( t ) . (69)The explicit solution of equation (64) in this case is: ξ N ( t ) = [ b N (0)] k N [ b N (0)] t + 1 . (70)Therefore, the equilibrium solution of equation (64) is:lim t →∞ ξ N ( t ) = ξ ∗ N = 0 . We obtained B N − and LD N − group concentration η N − ( t ) and the equilibrium value η ∗ N − , using themass conservation law equation (52). η N − ( t ) = [ b N − (0)] − [ b N (0)] + ξ N ( t ) . (71) η ∗ N − = [ b N (0)] − [ b N (0)] + ξ ∗ N . (72) A.3 LU reaction group Next, we considered LU reaction group. LU reaction group consisted of monomer X , , X , , · · · , X i,i · · · complexes, and X N,N . dξ ( t ) dt is obtained as the sum of the upper ODEs. The reaction terms definingwhich monomers or complexes of group LU lead to the formation of other monomers or complexes ofthis group are canceled. dξ ( t ) dt = − k ξ ( t ) + l { [ b (0)] − ξ ( t ) } = − k ξ ( t ) − l ξ ( t ) + l [ b (0)] . (73)The discriminant of dξ ( t ) dt is: D = l + 8 k l [ b (0)] . (74)If D is greater than 0, i.e., l >
0, there are two equilibrium value solutions of dξ ( t ) dt = 0: ξ ∗ ± = − l ± (cid:112) l + 8 k l [ b (0)]4 k . (75)We considered ξ ( t ) as the total monomer and complex concentrations in group LU : ξ ( t ) = ξ ∗ − C ξ ∗ − e − β t − C e − β t , (76)here C = [ b (0)]+ ξ ∗ − [ b (0)]+ ξ ∗ , β = k ( ξ ∗ − ξ ∗ − ). When time t tends to + ∞ , e − β t → +0:lim t →∞ ξ ( t ) = ξ ∗ . The discriminant D = 0 when l = 0. In this case, we obtained a steady solution for the equation (3.1), ξ ∗ = 0 and equation dξ ( t ) dt = − k ξ ( t ) from (3.1). The total monomer and complex concentrations ingroup LU were: ξ ( t ) = [ b (0)] k [ b (0)] t + 1 , and ξ ( t ) decreases with O ( t − ) more slowly than in the case D > O ( e − β t ). When time t tendsto + ∞ , the numerator is constant and denominator tends to + ∞ .lim t →∞ ξ ( t ) = ξ ∗ = 0 . We obtain B and LD group concentration η ( t ), and the equilibrium value η ∗ , through the massconservation law equation (52): η ( t ) = [ b (0)] − [ b (0)] + ξ ( t ) . (77) η ∗ = [ b (0)] − [ b (0)] + ξ ∗ . (78)The LU reaction group concentrations are summarized as follows: ξ N ( t ) = (cid:40) ξ ∗ N + − C N ξ ∗ n − e − βnt − C n e − βnt , ( l N > , m = N ) , [ b N (0)] k N [ b N (0)] t +1 , ( l N = 0 and [ b N − (0)] = [ b N (0)] , m = N ) . (79) ξ m ( t ) = (cid:40) ξ ∗ m + − C m ξ ∗ m − e − βmt − C m e − βmt , ( l m > , m : 1 ↓ N − , [ b m (0)] k m [ b m (0)] t +1 , ( l m = 0 and [ b m − (0)] = [ b m (0)] , m : 1 ↓ N − . (80) ξ ( t ) = ξ ∗ − C ξ ∗ − e − β t − C e − β t , ( l > , m = 1) , [ b (0)] k [ b (0)] t +1 , ( l = 0 , m = 1) . (81)LU group concentrations, ξ N , ..., ξ m , ..., ξ , are independent of each other and have strict solutions, whichconverge to equilibrium solutions when time t tends to + ∞ .The total concentration of B m and LD m groups η m is derived from the equation (7), as follows: (cid:40) η N ( t ) = [ b N (0)] , ( m = N ) ,η m ( t ) = [ b m (0)] − [ b m +1 (0)] + ξ m +1 ( t ) , ( m = 1 , · · · , N − . (82)Here, we showed that the total concentration of LU m group, or total concentrations of LD m and B m groups, can be solved. They converge to the equilibrium solutions when time t tends to + ∞ .We regarded the reaction group concentrations as a priori upper boundaries of the concentration ofeach molecular in the group. Afterward, we demonstrated that each molecule in a group converges to astable solution. eferences [1] H Sato, T Takino, Y Okada, J Cao, A Shinagawa, E Yamamoto, and M Seiki,A matrix metallo-proteinase expressed on the surface of invasive tumor cells, Nature, (1994), 61-65.[2] ED Karagiannis, and AS Popel, A theoretical model of type I collagen proteolysis by matrixmetalloproteinase (MMP) 2 and membrane type 1 MMP in the presence of tissue inhibitor ofmetalloproteinase 2, J Biol Chem, (2004), 39105-39114.[3] D Hoshino, N Koshikawa, T Suzuki, K Ichikawa, V Quaranta, A Weaver, and M Seiki, Estab-lishment and validation of computational model for MT1-MMP dependent ECM degradation andintervention strategies, PLoS Comput Biol, (2012), 33[5] S. Kawasaki, D. Minerva, K. Itano and T. Suzuki, Finding Solvable Units of Variables in NonlinearODEs of ECM Degradation Pathway Network, Comput. Math. Methods Med., (2017), 1-15.[6] K. Itano and T. Suzuki, Mathematical Modelling and Analysis of the Enhanced ECM DegradationPathway Network with N Monomers., JSIAM J., (2016), 44-83.[7] G Shinar and M Feinberg, Structural Sources of Robustness in Biochemical Reaction Networks, (2010), 1389-1391.[8] A Mochizuki and B Fiedler, Sensitivity of chemical reaction networks: A structural approach. 1.Examples and the carbon metabolic networks, J Theor Biol, (2015), 189-202.[9] M Eitoku, L Satoa , T Senda, and M Horikoshi, Histone chaperones: 30 years from isolation toelucidation of the mechanisms of nucleosome assembly and disassembly, Cell Mol Life Sci, (2008),414-444.[10] T Nogi, N Yasui, E Mihara, Y Matsunaga, M Noda, N Yamashita, T Toyofuku, S Uchiyama, YGoshima, A Kumanogoh, and J Takagi, Structural basis for semaphorin signalling through theplexin receptor, Nature, (2010), 1123-1127.[11] NN Noda, Y Fujioka, T Hanada, Y Ohsumi, and F Inagaki, Structural basis for semaphorinsignalling through the plexin receptor, Nature, (2010), 1123-1127.[12] DJ Klionsky and BA Schulman, Dynamic regulation of macroautophagy by distinctive ubiquitin-like proteins, Nat Struct Mol Biol, (2014), 336-345.[13] M Sakoh-Nakatogawa, H Kirisako, H Nakatogawa, and Y Ohsumi, Localization of Atg3 toautophagy-related membranes and its enhancement by the Atg8-family interacting motif to pro-mote expansion of the membranes, FEBS Lett, (2014), 336-345.[14] H Nakatogawa, Two ubiquitin-like conjugation systems that mediate membrane formation duringautophagy, Essays Biochem,55