# Finding Solvable Units of Variables in Nonlinear ODEs of ECM Degradation Pathway Network

Authors: Shuji Kawasaki; Dhisa Minerva; Keiko Itano; Takashi Suzuki
Journal: Computational and Mathematical Methods in Medicine
 (2017)

---

## Abstract

We consider ordinary differential equation (ODE) model for a pathway network that arises in extracellular matrix (ECM) degradation. For solving the ODEs, we propose applying the mass conservation law (MCL), together with a stoichiometry calleddoubling rule, to them. Then it leads to extracting new units of variables in the ODEs that can be solved explicitly, at least in principle. The simulation results for the ODE solutions show that the numerical solutions are indeed in good accord with theoretical solutions and satisfy the MALs.


---

## Body

## 1. Introduction

Differential equations are a useful method for modeling dynamics of reaction pathways in cells. They can be used to formulate biochemical kinetics, that is, the interactive dynamics of systems consisting of proteins, enzymes, products, and other components in terms of their transitive concentrations.The biochemical kinetics of a system describes the reactions through the mass action law (MAL); this results in a system of first-order nonlinear ordinary differential equations (ODEs). In general, the nonlinear ODEs cannot be solved explicitly, and so they are usually studied by numerical simulations.For the modeling of such a system, we propose an idea of using the mass conservation laws (MCLs) together with the MALs. That is, the variables in the ODEs obtained by the MALs are grouped into new units of variables to constitute new MALs, according to local balancing relations of inflows and outflows. It is important here that these new MCLs hold for allt≥0. As a result, the nonlinear ODEs presented by the new units of variables turns out to be acompletely integrable system; that is, the ODEs can be solved explicitly, at least in principle.We apply this idea to analyze the kinetics of the molecule concentration dynamics in cancer cell invasion to the ECM by a matrix metalloproteinase called MT1-MMP. If we can find appropriate MCLs for balancing relations in the kinetics, then such relations turn out to exhibit linear relations between ODE variables that are valid for allt≥0. Also, the nonlinear ODEs with the new unit of variables are completely solvable explicitly. That is, these nonlinear ODEs of new unit of variables can be solvedas a group. These groups may often correspond to network motifs, that is, local functions; they suggest a bundle of meaningful components in a complex pathway network (PWN). As the simulation results show, the numerical solutions are in good accord with the theoretical solutions, indicating our modeling by the new unit of variables is right.We thus clarify how to look at the PWN model of the ECM degradation with the unit variables. In fact, we indicate that the relevant system of the PWN can be grouped into several units of original variables and, with the unit variables, the system of ODEs can be solved explicitly. Taking linear combination of ODEs for the grouping has existed so far, but changing original ODEs into new, solvable ones by the linear combination has not been considered so far. This is a special stoichiometry. In addition, in order to enjoy such grouping method, certain reformulation includingdoubling rule is necessary as in the following, which has not been considered properly so far.The organization of the paper is as follows. Section2 discusses four elementary reaction processes that are modeled by second-order nonlinear ODEs that can be solved explicitly. In Section 3, the ECM degradation mechanism is introduced, and we present one of the key concepts behind the choice of a good model for ODE kinetics, thedoubling rule in stoichiometry. In Section 4, as an application of elementary reaction processes, we consider solving the ODE system that arises from the PWN of molecule concentration dynamics in ECM degradation. We list the tables cited in Section 4 and present the explicit solutions to the ODE system in Tables 1–7 and Appendix A, respectively.Table 1
List of associations witha-bB.
Association witha-bB
Assoc. const.
Dissoc. const.
X
1
(
a
)
+
X
2
b
⟶
X
4
a
b
k
1
l
1
=
0
+
X
5
b
c
⟶
X
7
a
b
c
k
1
+
X
8
(
b
c
c
)
⟶
X
9
(
a
b
c
c
)
k
1
+
X
10
(
b
c
c
b
)
⟶
X
11
(
a
b
c
c
b
)
2
k
1
+
X
11
(
a
b
c
c
b
)
⟶
X
12
(
a
b
c
c
b
a
)
k
1Table 2
List of associations/dissociations withBb-cC.
Association/dissociation withBb-cC
Assoc. const.
Dissoc. const.
X
2
(
b
)
+X3(c)⟷X5(bc)
k
2
l
2
+X6(cc)⟷X8(bcc)
2
k
2
l
2
+X8(bcc)⟷X10(bccb)
k
2
2
l
2
+X9(abcc)⟷X11(abccb)
k
2
l
2
X
4
(
a
b
)
+X3(c)⟷X7(abc)
k
2
l
2
+X6(cc)⟷X9(abcc)
2
k
2
l
2
+X8(bcc)⟷X11(abccb)
k
2
l
2
+X9(abcc)⟷X12(abccba)
k
2
2
l
2Table 3
List of associations/dissociations withCc-cC.
Association/dissociation withCc-cC
Assoc. const.
Dissoc. const.
X
3
(
c
)
+
X
3
(
c
)
⟷
X
6
(
c
c
)
2
k
3
2
l
3
+
X
5
(
b
c
)
⟷
X
8
(
b
c
c
)
k
3
l
3
+
X
7
(
a
b
c
)
⟷
X
9
(
a
b
c
c
)
k
3
l
3
X
5
(
b
c
)
+
X
3
(
c
)
⟷
X
8
(
b
c
c
)
k
3
l
3
+
X
5
(
b
c
)
⟷
X
10
(
b
c
c
b
)
2
k
3
2
l
3
+
X
7
(
a
b
c
)
⟷
X
11
(
a
b
c
c
b
)
k
3
l
3
X
7
(
a
b
c
)
+
X
3
(
c
)
⟷
X
9
(
a
b
c
c
)
k
3
l
3
+
X
5
(
b
c
)
⟷
X
11
(
a
b
c
c
b
)
k
3
l
3
+
X
7
(
a
b
c
)
⟷
X
12
(
a
b
c
c
b
a
)
2
k
3
2
l
3Table 4
Mass conservation laws ofa, b, and c.
Items
MCL
a
0
d
X
1
+
X
4
+
X
7
+
X
9
+
X
11
+
2
X
12
d
t
≡
0
b
0
d
X
2
+
X
4
+
X
5
+
X
7
+
X
8
+
X
9
+
2
X
10
+
2
X
11
+
2
X
12
d
t
≡
0
c
0
d
X
3
+
X
5
+
2
X
6
+
X
7
+
2
X
8
+
2
X
9
+
2
X
10
+
2
X
11
+
2
X
12
d
t
≡
0Table 5
Reaction laws of(a,b), (b,c), and (c,c).
a
,
b: X1+X2+X5+X8+2X10+X11⇌l1≡0k1X4+X7+X9+X11+2X12⟺
E
1
d
X
1
d
t
=
-
k
1
X
1
X
2
+
X
5
+
X
8
+
2
X
10
+
X
11
E
2581011
d
X
2
+
X
5
+
X
8
+
2
X
10
+
X
11
d
t
=
-
k
1
X
1
X
2
+
X
5
+
X
8
+
2
X
10
+
X
11
.
b
,
c: X2+X4+X3+2X6+X8+X9⇌l2k2X5+X7+X8+X9+2X10+2X11+2X12⟺
E
24
d
X
2
+
X
4
d
t
=
-
k
2
X
2
+
X
4
X
3
+
2
X
6
+
X
8
+
X
9
+
l
2
X
5
+
X
7
+
X
8
+
X
9
+
2
X
10
+
2
X
11
+
2
X
12
E
3689
d
X
3
+
2
X
6
+
X
8
+
X
9
d
t
=
-
k
2
X
2
+
X
4
X
3
+
2
X
6
+
X
8
+
X
9
+
l
2
X
5
+
X
7
+
X
8
+
X
9
+
2
X
10
+
2
X
11
+
2
X
12
.
(
c
,
c
): X3+X5+X7+X3+X5+X7⇌l3k32X6+X8+X9+X8+2X10+X11+X9+X11+2X12⟺
(
E
357
)
d
X
3
+
X
5
+
X
7
d
t
=
-
2
k
3
X
3
+
X
5
+
X
7
2
+
2
l
3
X
6
+
X
8
+
X
9
+
X
10
+
X
11
+
X
12Table 6
Processes of obtaining solutionsX1–X12.
Variable
Method
X
1
Single Riccati equation
ξ
2581011
A companion ofX1
ξ
24
Single Riccati equation
ξ
3689
A companion ofξ24
X
2
Method of variation of coefficient
X
4
X
4
=
ξ
24
-
X
2
ξ
357
Single Riccati equation
X
3
,
X
5
Method of variation of coefficient
X
7
X
7
=
ξ
357
-
(
X
3
+
X
5
)
X
6, X8, X10
Method of variation of coefficient
X
9, X11, X12
by (41)Table 7
Initial values and reaction constants used in the simulation.
Variables
Values
a
0
2.0
×
1
0
-
6 [M]
b
0
1.0
×
1
0
-
6 [M]
c
0
1.0
×
1
0
-
6 [M]
k
1
2.1
×
1
0
7 [M-1s-1]
l
1
0 [M-1s-1]
k
2
2.74
×
1
0
6 [M-1s-1]
l
2
2.0
×
1
0
-
4 [M-1s-1]
k
3
2.0
×
1
0
6 [M-1s-1]
l
3
1.0
×
1
0
-
2 [M-1s-1]


## 2. Elementary Reaction Processes That Reduce to Solvable Second-Order Riccati Equations

In this section, we will show how the ODEs for the elementary reaction processes can be solved by reduction to the well-known Riccati equations of a new unit of variables.In PWN, in addition to the simple dimerization of monomers, there are polymerizations of higher complexes. For those complexes that are modified from relevant monomers, we will assume that the association and dissociation constants are the same as those for the monomer associations or dimer dissociations, respectively.For simplicity, we will show the four basic forms of elementary reactions as in the following. The association dissociation constants are here denoted byk and l, respectively.(i)
Association of two different molecules/dissociation of their product:(1)A+B⇌lkAB.(ii)
Association of the same molecules/dissociation of their product:(2)A+A⇌lkAA.(iii)
Association with a modified molecule/dissociation of their product:(3)A+AB⇌lkABAorAAB.(iv)
Association with a homodimer/dissociation of their product:(4)A+BB⇌lkABB.Note that AB+B⇌lkABB is the same as reaction (iii).We now show that the ODEs for the MALs for reactions (i)–(iv) can be solved explicitly. First, the MALs for reaction (i) are(5)dAdt=-kAB+lP,dBdt=-kAB+lP,dPdt=kAB-lP.The resulting system of ODEs can be solved explicitly: (5) implies the MCLs(6a)dA+Pdt=dB+Pdt=0⟺A+P≡IA,B+P≡IB,(6b)dA-Bdt=0⟺A-B=IA-IBfor some positive constants IA and IB, which are typically the initial values of [A]+[P] and [B]+[P], respectively. Then, substituting (7)A=IA-P,B=IB-Pinto the third equation in (5), we obtain an equation in only the variable [P]: (8)dPdt=kIA-PIB-P-lP=kP2-kIA+IB+lP+kIAIB≡kP-z1P-z2,where the last expression is the factorization of the quadratic polynomial in the second line. If z1≠z2, this leads to (9)1P-z2-1P-z1dPdt=-kz1-z2,or P-z1/P-z2=Ce-k(z1-z2)t, and hence the explicit solution (10)Pt=z1-Cz2e-kz1-z2t1-Ce-kz1-z2t=z2+z1-z21-Ce-kz1-z2t,C=P0-z1P0-z2for t≥0. Here, z1,z2 are (11)z1,z2=kIA+IB+l±kIA+IB+l2-4k2IAIB2k.If z1 and z2 are distinct and z1-z2>0 so that limt→∞[P]t=z1 in (10), then we find that [P]t is monotonically decreasing and nonnegative. Solutions to [A] and [B] can be found in a similar manner by using (7). Although the situation is slightly different, a similar calculation was presented in ([1], Section  6.5).Second, for reaction (ii), the MAL yields kinetic equations of the form(12)dAdt=-2kA2+2lP,dPdt=kA2-lP.In the first equation of (12), the coefficient of 2 stems from the fact that, in the reaction A+A→P, two molecules of A are consumed; similarly, in the dissociation, two molecules of A are produced. This is called adouble rate of consumption/reproduction. The ODEs (12) imply the MCL (13)dA+2Pdt=0⟺A+2P≡IAfor a positive constant IA. Then, substituting 2[P]=IA-[A] into the first equation of (12), we have (14)dAdt=-2kA2-lA+lIA≡-2kA-ζ1A-ζ2,At=ζ1-Cζ2e-2kζ1-ζ2t1-Ce-2kζ1-ζ2t=ζ2+ζ1-ζ21-Ce-2kζ1-ζ2t,C=A0-ζ1A0-ζ2for t≥0. Here ζ1,ζ2 are (15)ζ1,ζ2=-l±l2+8klIA4k.Since ζ1-ζ2>0, we have limt→∞[A]t=ζ1. As in case (i), the solution is monotonically decreasing and nonnegative. The solution to [P] can also be obtained by the MCL (13).For reaction (iii), we have the MALs and MCL as shown in the following:(16)dAdt=-kAAB+lP,dABdt=-kAAB+lP,dPdt=kAAB-lP,(17a)dA+Pdt=0,(17b)dA-ABdt=0.The ODEs (16) can be solved in a manner similar to that used for reactions (i) and (ii). In fact, the association of A-BA reduces to reaction (i), and the association of A-AB reduces to reaction (ii).Finally, for reaction (iv), we have the following MALs and the MCL:17dAdt=-2kABB+lABB,dBBdt=-2kABB+lABB,dABBdt=2kABB-lABB,(18)dA-BBdt=0,dA+ABBdt=0.The coefficient of 2 in 17 arises because a molecule of A may bind to either of the two B molecules. The solution to these ODEs 17 is the same as that for reaction (i) but with k replaced by 2k. This is called adouble chance of association.These elementary reaction processes (i)–(iv) will be used below to solve the ODE system for the kinetics of ECM degradation. The given ODE system is grouped into several subpathways, and in the reformulation based on these groups, we will find that the new variables reduce the equationsas groups to those of the elementary reaction processes.


## 3. Application to the Model of ECM Degradation by MT1-MMP

In this section, we show how the elementary reaction processes in the previous section can be applied to a problem in cell biology, that of extracellular matrix (ECM) degradation by membrane type 1 matrix metalloproteinase (MT1-MMP), which is a step in the progression of cancer.


### 3.1. ECM Degradation Mechanism

As is well known, in order for the proliferated cancer cells to metastasize from a primary lesion, they first invade and degrade the ECM and then migrate. The cancer cells secrete MMP, and its dimer MMP2 is then activated by MT1-MMP; the MMP2 then degrades the ECM.After the ECM is degraded, other enzymes, including MT1-MMP, begin degrading the interstitium beyond the ECM. Therefore, mathematical or biological elucidation of the activation process of MMP2 is important for finding a clinical cure or developing drugs. Below, we will denote MMP2, TIMP2, and MT1-MMP bya, b, and c, respectively.Sato et al. [2] have revealed experimentally the mechanism of the activation of MMP2. The steps of the activation scenario of MMP2, as a cell biological model, proceed as follows (see Figure 1):(0)
Althoughc (MT1-MMP) is an ECM degradation enzyme as well, it is usually inactivated by b (TIMP2) immediately after vesicle transportation to the cell membranes. However, c activates a (MMP2), which is another ECM degradation enzyme, as follows.(1)
Thec molecules that penetrate the cell membrane form homodimers cc on the membrane.(2)
On one of thec in the dimer -cc-, the heterodimer ab (pro-MMP2, TIMP2) is coupled to produce the tetramer abcc; alternatively, a- may be coupled with the heterotrimer -bcc to produce the tetramer. Thus, a is associated with cc via b.(3)
One of thec in abcc that is not coupled with ab cuts the connection between a and b.(4)
a that was cut and released then becomes the activated one.This process is thus the associations/dissociations of the three molecules a, b, and c, and abcc is the key molecule in the activation of a. In order to produce certain amount of activated a (MMP2), there must be sufficient b. However, if there is too much b, then the remaining vacant site of abcc- tends to associate with -b, which results in the fact that abccb or abccba will be produced preferentially to abcc, and an insignificant amount of a will be activated. See Hoshino et al. [3] and then Watanabe et al. [4].Figure 1
Mechanism of activation ofa (MMP2) byb (TIMP2) andc (MT1-MMP).The molecular binding rules that follow from this scenario are thata and c do not couple with each other, but with b, and c may form the homodimer cc (see Figure 2). Although it is difficult to deny the possibility of other binding patterns, to the best of our knowledge, there is no evidence of such patterns in the literature. Therefore, we will assume the following binding patterns:(1)
a has a site of binding only with b.(2)
b has a site of binding with a and another site of binding with c.(3)
c has a site of binding with b and another site of binding with c itself.(4)
There are no other binding patterns.Figure 2
Binding sites ofa (MMP2),b (TIMP2), andc (MT1-MMP).By inspection, we see that there appear to be nine complexes (dimers to hexamer) obtainable from the monomersa, b, and c. The PWN of association or dissociation of these complexes are shown in Figure 3.Figure 3
Pathway network of association/dissociation pairs of complexes ofa (MMP2),b (TIMP2), andc (MT1-MMP). Here, arrows indicating only association are depicted.Thus, a quantitative model of ECM degradation can be based on the kinetics of associations/dissociations of the three proteins (a: MMP2; b: TIMP2; and c: MT1-MMP) and their nine complexes, for a total of twelve compounds. The PWNs of the twelve molecules are depicted in Figure 3. A quantitative model was first considered by Karagiannis and Popel [5], in which the general interactive behavior of the molecules a, b, and c was investigated through simulation in the context of type-I collagen proteolysis. Further, their behavior in more realistic situations was considered by Hoshino et al. [3] and then by Watanabe et al. [4] as the mathematical interpretation of the cell biological model in [2].In the quantitative models in [3, 4], the detailed cell biological behaviors of these enzymes (a, b, and c) were studied. In particular, they found that MT1-MMP has two fluorescence recovery time constants, 29 [s] and 259 [s], and, in FRAP experiments, they determined that the recovery is not due to lateral diffusions but to the insertion of a new membrane. As a result, they verified by simulation that the inactivation of MT1-MMP by TIMP2 proceeds very quickly (halving time is 4.5 [s]); therefore, a very rapid turnover of MT1-MMP must be occurring at the location of the ECM degradation.We consider a mathematical treatment of the ODE kinetics with a more complete theoretical analysis of such quantitative studies; seeE1–E12.Now, according to the binding rules discussed above, we have the following three associations or dissociations in the PWN:(I)
Association ofa-bB (association constant k1): (19)a-ba+bB⇌l1≡0k1abB.The molecules a and b do not dissociate after they are once associated. Thus we set l1≡0.(II)
Association/dissociation ofBb-cC (association constant k2; dissociation constant l2): (20)b-cbB+cC⇌l2k2BbcC.(III)
Association/dissociation ofCc-cC (association constant k3; dissociation constant l3): (21)c-ccC+cC⇌l3k3CccC.For reactions (I)–(III), there are twelve molecules in the PWN in Figure3: a, b, c, ab, bc, cc, abc, bcc, abcc, bccb, abccb, and abccba. We denote their concentrations as X1=[a], X2=[b], X3=[c], and so on. The three association/dissociation groups are presented in Tables 1, 2, and 3. We make the following basic assumptions:(A1)
For the association/dissociation of the modified moleculesAa and bB, we use the same association/dissociation constants as those used for the elementary reaction a-b. The same is assumed for Bb-cC and Cc-cC.(A2)
The initial concentrations forX4,…,X12 are all 0: Xi(0)=0, i=4,…,12.In order to obtain the correct ODEs, it is necessary to extract the appropriatedoubling rules for the reaction coefficients in the PWN. We will explain it in the next subsection.


### 3.2. Pathway Network Dynamics and Doubling Rules

From Figure3 and Tables 1–3, we can write down the following twelve second-order nonlinear ODEs:E1dX1dt=-k1X1X2+X5+X8+2X10+X11.E2dX2dt=-k1X1X2-k2X2X3+2X6+X8+X9+l2X5+X8+2X10+X11.E3dX3dt=-k2X3X2+X4-2k3X3X3+X5+X7+l2X5+X7+l32X6+X8+X9.E4dX4dt=k1X1X2-k2X4X3+2X6+X8+X9+l2X7+X9+X11+2X12.E5dX5dt=-k1X1X5+k2X2X3-2k3X5X3+X5+X7-l2X5+l3X8+2X10+X11.E6dX6dt=k3X32-2k2X6X2+X4+l2X8+X9-l3X6.E7dX7dt=k1X1X5+k2X3X4-2k3X7X3+X5+X7-l2X7+l3X9+X11+2X12.E8dX8dt=-k1X1X8+2k2X2X6-k2X8X2+X4+2k3X3X5+l2-X8+2X10+X11-l3X8.E9dX9dt=k1X1X8-k2X9X2+X4+2k2X4X6+2k3X3X7+l2-X9+X11+2X12-l3X9.E10dX10dt=-2k1X1X10+k2X2X8+k3X52-2l2X10-l3X10.E11dX11dt=k1X12X10-X11+k2X2X9+k2X4X8+2k3X5X7-2l2X11-l3X11.E12dX12dt=k1X1X11+k2X4X9+k3X72-2l2X12-l3X12.Note that ODEsE1–E12 incorporate doubling rules. The coefficient of 2 on(1)
k3X32 in E3,(2)
k3X52 in E5,(3)
k3X72 in E7corresponds to the double rate of consumption/reproduction in reaction (ii). See Figure 4. The coefficient of 2 on (1)
k3X3(X5+X7) in E3,(2)
k3X5(X3+X7) in E5,(3)
k3X7(X3+X5) in E7,(4)
k3X3X5 in E8,(5)
k3X3X7 in E9,(6)
k3X5X7 in E11,which corresponds to Cc-cC association or dissociation reactions in Table 3, will be explained in Section 4.3. The remaining coefficients of 2 in the ODEs are due to the double chance of association/dissociation in reaction (iv), as discussed in Section 2.Figure 4
Doubling rules in the association/dissociation ofb- and -cC.What is the rationale for the above doubling rules? It is the ingeniously well-formed MCLs and reaction laws, which are shown in the next section.


## 3.1. ECM Degradation Mechanism

As is well known, in order for the proliferated cancer cells to metastasize from a primary lesion, they first invade and degrade the ECM and then migrate. The cancer cells secrete MMP, and its dimer MMP2 is then activated by MT1-MMP; the MMP2 then degrades the ECM.After the ECM is degraded, other enzymes, including MT1-MMP, begin degrading the interstitium beyond the ECM. Therefore, mathematical or biological elucidation of the activation process of MMP2 is important for finding a clinical cure or developing drugs. Below, we will denote MMP2, TIMP2, and MT1-MMP bya, b, and c, respectively.Sato et al. [2] have revealed experimentally the mechanism of the activation of MMP2. The steps of the activation scenario of MMP2, as a cell biological model, proceed as follows (see Figure 1):(0)
Althoughc (MT1-MMP) is an ECM degradation enzyme as well, it is usually inactivated by b (TIMP2) immediately after vesicle transportation to the cell membranes. However, c activates a (MMP2), which is another ECM degradation enzyme, as follows.(1)
Thec molecules that penetrate the cell membrane form homodimers cc on the membrane.(2)
On one of thec in the dimer -cc-, the heterodimer ab (pro-MMP2, TIMP2) is coupled to produce the tetramer abcc; alternatively, a- may be coupled with the heterotrimer -bcc to produce the tetramer. Thus, a is associated with cc via b.(3)
One of thec in abcc that is not coupled with ab cuts the connection between a and b.(4)
a that was cut and released then becomes the activated one.This process is thus the associations/dissociations of the three molecules a, b, and c, and abcc is the key molecule in the activation of a. In order to produce certain amount of activated a (MMP2), there must be sufficient b. However, if there is too much b, then the remaining vacant site of abcc- tends to associate with -b, which results in the fact that abccb or abccba will be produced preferentially to abcc, and an insignificant amount of a will be activated. See Hoshino et al. [3] and then Watanabe et al. [4].Figure 1
Mechanism of activation ofa (MMP2) byb (TIMP2) andc (MT1-MMP).The molecular binding rules that follow from this scenario are thata and c do not couple with each other, but with b, and c may form the homodimer cc (see Figure 2). Although it is difficult to deny the possibility of other binding patterns, to the best of our knowledge, there is no evidence of such patterns in the literature. Therefore, we will assume the following binding patterns:(1)
a has a site of binding only with b.(2)
b has a site of binding with a and another site of binding with c.(3)
c has a site of binding with b and another site of binding with c itself.(4)
There are no other binding patterns.Figure 2
Binding sites ofa (MMP2),b (TIMP2), andc (MT1-MMP).By inspection, we see that there appear to be nine complexes (dimers to hexamer) obtainable from the monomersa, b, and c. The PWN of association or dissociation of these complexes are shown in Figure 3.Figure 3
Pathway network of association/dissociation pairs of complexes ofa (MMP2),b (TIMP2), andc (MT1-MMP). Here, arrows indicating only association are depicted.Thus, a quantitative model of ECM degradation can be based on the kinetics of associations/dissociations of the three proteins (a: MMP2; b: TIMP2; and c: MT1-MMP) and their nine complexes, for a total of twelve compounds. The PWNs of the twelve molecules are depicted in Figure 3. A quantitative model was first considered by Karagiannis and Popel [5], in which the general interactive behavior of the molecules a, b, and c was investigated through simulation in the context of type-I collagen proteolysis. Further, their behavior in more realistic situations was considered by Hoshino et al. [3] and then by Watanabe et al. [4] as the mathematical interpretation of the cell biological model in [2].In the quantitative models in [3, 4], the detailed cell biological behaviors of these enzymes (a, b, and c) were studied. In particular, they found that MT1-MMP has two fluorescence recovery time constants, 29 [s] and 259 [s], and, in FRAP experiments, they determined that the recovery is not due to lateral diffusions but to the insertion of a new membrane. As a result, they verified by simulation that the inactivation of MT1-MMP by TIMP2 proceeds very quickly (halving time is 4.5 [s]); therefore, a very rapid turnover of MT1-MMP must be occurring at the location of the ECM degradation.We consider a mathematical treatment of the ODE kinetics with a more complete theoretical analysis of such quantitative studies; seeE1–E12.Now, according to the binding rules discussed above, we have the following three associations or dissociations in the PWN:(I)
Association ofa-bB (association constant k1): (19)a-ba+bB⇌l1≡0k1abB.The molecules a and b do not dissociate after they are once associated. Thus we set l1≡0.(II)
Association/dissociation ofBb-cC (association constant k2; dissociation constant l2): (20)b-cbB+cC⇌l2k2BbcC.(III)
Association/dissociation ofCc-cC (association constant k3; dissociation constant l3): (21)c-ccC+cC⇌l3k3CccC.For reactions (I)–(III), there are twelve molecules in the PWN in Figure3: a, b, c, ab, bc, cc, abc, bcc, abcc, bccb, abccb, and abccba. We denote their concentrations as X1=[a], X2=[b], X3=[c], and so on. The three association/dissociation groups are presented in Tables 1, 2, and 3. We make the following basic assumptions:(A1)
For the association/dissociation of the modified moleculesAa and bB, we use the same association/dissociation constants as those used for the elementary reaction a-b. The same is assumed for Bb-cC and Cc-cC.(A2)
The initial concentrations forX4,…,X12 are all 0: Xi(0)=0, i=4,…,12.In order to obtain the correct ODEs, it is necessary to extract the appropriatedoubling rules for the reaction coefficients in the PWN. We will explain it in the next subsection.


## 3.2. Pathway Network Dynamics and Doubling Rules

From Figure3 and Tables 1–3, we can write down the following twelve second-order nonlinear ODEs:E1dX1dt=-k1X1X2+X5+X8+2X10+X11.E2dX2dt=-k1X1X2-k2X2X3+2X6+X8+X9+l2X5+X8+2X10+X11.E3dX3dt=-k2X3X2+X4-2k3X3X3+X5+X7+l2X5+X7+l32X6+X8+X9.E4dX4dt=k1X1X2-k2X4X3+2X6+X8+X9+l2X7+X9+X11+2X12.E5dX5dt=-k1X1X5+k2X2X3-2k3X5X3+X5+X7-l2X5+l3X8+2X10+X11.E6dX6dt=k3X32-2k2X6X2+X4+l2X8+X9-l3X6.E7dX7dt=k1X1X5+k2X3X4-2k3X7X3+X5+X7-l2X7+l3X9+X11+2X12.E8dX8dt=-k1X1X8+2k2X2X6-k2X8X2+X4+2k3X3X5+l2-X8+2X10+X11-l3X8.E9dX9dt=k1X1X8-k2X9X2+X4+2k2X4X6+2k3X3X7+l2-X9+X11+2X12-l3X9.E10dX10dt=-2k1X1X10+k2X2X8+k3X52-2l2X10-l3X10.E11dX11dt=k1X12X10-X11+k2X2X9+k2X4X8+2k3X5X7-2l2X11-l3X11.E12dX12dt=k1X1X11+k2X4X9+k3X72-2l2X12-l3X12.Note that ODEsE1–E12 incorporate doubling rules. The coefficient of 2 on(1)
k3X32 in E3,(2)
k3X52 in E5,(3)
k3X72 in E7corresponds to the double rate of consumption/reproduction in reaction (ii). See Figure 4. The coefficient of 2 on (1)
k3X3(X5+X7) in E3,(2)
k3X5(X3+X7) in E5,(3)
k3X7(X3+X5) in E7,(4)
k3X3X5 in E8,(5)
k3X3X7 in E9,(6)
k3X5X7 in E11,which corresponds to Cc-cC association or dissociation reactions in Table 3, will be explained in Section 4.3. The remaining coefficients of 2 in the ODEs are due to the double chance of association/dissociation in reaction (iv), as discussed in Section 2.Figure 4
Doubling rules in the association/dissociation ofb- and -cC.What is the rationale for the above doubling rules? It is the ingeniously well-formed MCLs and reaction laws, which are shown in the next section.


## 4. Analyzing the PWN Dynamics

### 4.1. Kinetics of Those Reactions Containinga-bB

According to Table1, the reactions related to a-bB are summarized as (22)X1+X2+X5+X8+2X10+X11⇌l1≡0k1X4+X7+X9+X11+2X12.Looking atX1 and the products (the right hand side) of this reaction, we infer from the flux balance of consumption-production that we may have (23)dX1+X4+X7+X9+X11+2X12dt≡0.By taking summation E1 + E4 + E7 + E9 + E11 + 2 × E12 of the ODEs, we find that (23) is indeed true. Here, X12(abccba) has a coefficient of 2 because the decrease in X1 is equal to the increase in (X4+X7+X9+X11+2X12), since the complex abccba consumes two a molecules, unlike the others: ab,abc,abcc,abccb. Equation (23) is an MCL of form (6a) or (17a). Equation (23) may be considered as an MCL for a0, since (24)X1+X4+X7+X9+X11+2X12t≡const.≡X10=a0by assumption (A2); thus, a local balance holds: a0-X1(t)≡(X4+X7+X9+X11+2X12)(t) for all t≥0.From (22), we also have (25)d-X1+X2+X5+X8+2X10+X11dt≡0,which may be considered to represent an MCL for a0-b0, since (26)-X1+X2+X5+X8+2X10+X11t≡const.≡-X10+X20=-a0+b0.Here, X10 has a coefficient of 2 because the complex bccb has two b- sites for binding with a, which implies a double chance of association, as in reaction (iv). Equation (26) is an MCL of form (6b) or (17b).Note thatE1 and (25) suggest that X1 and (X2+X5+X8+2X10+X11) form a pair of a-bB kinetics, (E1) and (E2581011), as displayed in Table 5. The two ODEs can be solved explicitly using the method for reaction (i) in Section 2; details are shown in Appendix A.1.


### 4.2. Kinetics of Those Reactions ContainingBb-cC

Next, we consider the kinetics of associations/dissociations withBb-cC involved in the ODE system. As can be seen in Table 2, the reactions related to Bb-cC are summarized as (27)X2+X4+X3+2X6+X8+X9⇌l2k2X5+X8+2X10+X11+X7+X9+X11+2X12.In the manner similar to that used for the Aa-bB kinetics, we infer, from the flux balance of consumption-production, that we may have (28)dX2+X5+X8+2X10+X11+X4+X7+X9+X11+2X12dt=0;that is, (29)dX2+X4+X5+X7+X8+X9+2X10+2X11+2X12dt=0.By taking summation E1 + E4 + E7 + E9 + E11 + 2 × E12 of the ODEs, we find that (28) is indeed valid. Equation (28) is an MCL of form (6a) or (17a). Equstion (29) may be viewed as an MCL of b0, since (30)X2+X4+X5+X7+X8+X9+2X10+2X11+2X12t≡const.≡X20=b0.From (27), we also have (31)d-X2-X4+X3+2X6+X8+X9dt=0,which is confirmed by taking the summation −E2 − E4 + E3 + 2 × E6 + E8 + E9 of the ODEs. This is an MCL of form (6b) or (17b). Equation (31) represents an MCL of b0-c0, since (32)-X2-X4+X3+2X6+X8+X9t≡const.≡-X20+X30=-b0+c0.Here, X6 has a coefficient of 2 because it has two c- sites for binding with a or ab.The MCL (31) implies that ξ24(t)≜(X2+X4)(t) and ξ3689(t)≜(X3+2X6+X8+X9)(t) are governed by the same ODE, up to the initial values. From the ODEs E1–E12, it follows that the two variables form a pair of Bb-cC kinetics, (E24) and (E3689), as in Table 4. The variables ξ24(t) and ξ3689(t) can be solved explicitly; details are shown in Appendix A.2.Now, we find that each ofX2(t) and X4(t) can be solved explicitly, as well. This is because the right-hand side of E2 contains only known terms, except for X2: E2′dX2dt=-k1X1t+k2ξ24t+l2-k2b0-c0X2t+l2X1t-a0+b0,by (26) and (32). Hence X2(t) can be calculated by the method of variation of coefficients. See Appendix A.3. In this way, X4(t) can be obtained by X4(t)=ξ24(t)-X2(t).


### 4.3. Kinetics of Those Reactions ContainingCc-cC

We finally consider the kinetics of associations/dissociations ofCc-cC involved in the ODE system. As shown in Table 2, the reactions related to Cc-cC can be summarized as (33)X3+X5+X7+X3+X5+X7⇌l3k32X6+X8+X9+X8+2X10+X11+X9+X11+2X12.In a manner similar to that used previously, we infer from the flux balance of consumption-production that(34)dX3+2X6+X8+X9+X5+X8+2X10+X11+X7+X9+X11+2X12dt=0;that is, (35)dX3+X5+2X6+X7+2X8+2X9+2X10+2X11+2X12dt=0.By taking the summation E3 + E5 + 2 × E6 + E7 + 2 × E8 + 2 × E9 + 2 × E10 + 2 × E11 + 2 × E12 of the ODEs, we find that (35) is indeed valid. Equation (34) is an MCL of form (6a) or (17a). Equation (35) may be viewed as an MCL of c0, since (36)X3+X5+2X6+X7+2X8+2X9+2X10+2X11+2X12t≡const.≡X30=c0.From E3, E5, and E7, we have the ODE (E357), which is displayed in Table 5, and in which we have used (30) and (32). The coefficient of 2 on k3X3(X3+X5+X7)2 in (E357) is a result of the summation of the ODEs E3, E5, and E7. But then, where do the coefficients of 2 on k3(X3+X5+X7) in E3, on k3X5(X3+X5+X7) in E5, and on k3X7(X3+X5+X7) in E7 come from? They stem from the following stoichiometry. For those reactions in Table 3 that correspond to the reaction (ii) in Section 2, (37)X3c+X3c⟷X6cc,X5bc+X5bc⟷X10bccb,X7abc+X7abc⟷X12abccba,the reason is the same as in (12); that is, it is caused by the double rate of consumption/reproduction. Hence, we must count these three reactions with the association/dissociation coefficients on 2k3 and 2l3 in E3, E5, and E7. For the reaction X3(c)+X5(bc)↔X8(bcc) corresponding to reaction (iii), both (38)X3c+X5bc⟷X8bcc,X5bc+X3c⟷X8bccmust be counted. The first reaction in (38) is the one contained in E3, and the second reaction is in E5. In the first reaction, thesubject is X3c, and it associates with theobjectX5(bc). We will name this molecule c=c1. In the second reaction, thesubject is X5bc, and it associates with theobjectX3(c). We will name this molecule c=c2. In other words, a molecule of c (=object) is associated with a molecule of bc (=subject). In general, c=c1 and c2 are different molecules.The first reaction must be counted as a reaction ofX3c. Also, seen as a reaction of X5bc, the second reaction must also be counted. Both reactions take place at the same time. The association/dissociation coefficients are k3 and l3, respectively, for both reactions.However, now consider, for example, howX3(c) is consumed doubly by the amount of c1 and c2. In Table 3, other combinations of this kind are (39)X5bc+X7abc⟷X11abccb,X7abc+X5bc⟷X11abccb,X3c+X7abc⟷X9abcc,X7abc+X3c⟷X9abcc.The contribution of the reactions in (38) and (39) are therefore doubled. Combining the doubling of the reactions in (38), we obtain the third reaction in Table 5, which causes the term -2k3(X3+X5+X7)2 in E357 and the ODEs E3, E5, and E7.Finally, in consideringX3(c), the dissociation term l3(X8+X9) need not be doubled, since only one molecule of c is produced from a molecule of bcc: bc-c; also, this dissociation is not involved in E5.By the above arguments for the doubling rules, we have obtained ingeniously well-formed MCLs and reaction laws, as shown in Table4, respectively. This may suggest that our doubling rules are justified.Also, the above argument ofCc-cC kinetics may imply the following: suppose that, for association of molecules, p1,p2,…,pn, say, some combinations of coupling pi-pj (1≤i,j≤n), are possible and others not. Then, it does not hold that (40)dX1+⋯+Xndt=-2kX1+⋯+Xn2+dissociationtermin general. However, if all of pi has a common binding cite, that is, all pi are of the form -cC, then we do have (40). In the latter case, the doubling rule holds for (X1+⋯+Xn)as a unit, as if (p1,p2,…,pn) are a single molecule with binding cite -c, and the elementary doubling rule (12) holds with [A] in (12) replaced by (X1+⋯+Xn).


### 4.4. Solutions to Remaining Variables

Solutions to the remaining variables can also be obtained.X3 and X5 can be found by the method of variation of parameters, as in the case of X2. X7 can be found by X7(t)=ξ357(t)-(X3+X5)(t), where ξ357(t) is the solution to (X3+X5+X7)(t). By applying the solution ξ357 to E3, we obtain the solution to X3(t) as E3′dX3dt=-k2-l3ξ24+2k3-l2ξ357+l2+l3X3t+l2ξ357+l3ξ24-l3b0-c0by (30) and (32); note that ξ24 and ξ357 have been obtained already. Similarly, X6, X8, and X10 can also be found by the method of variation of coefficient.All ofX1 to X8 and X10 have thus been obtained in principle. The remaining ones, X9, X11, and X12, can be obtained by solving (41)X9t=ξ3689t-X3+2X6+X8t,X11t=X1t-X2+X5+X8+2X10t-a0+b0,X12t=12a0-X1+X4+X7+X9+X11t(see (A.8), (26), and (24), resp.). The processes for obtaining these solutions are listed in Table 6.


### 4.5. Biomedical Implication

As described in Section3.1, an appropriate amount of X2(b) (TIMP2) is necessary to obtain a substantial amount of the key molecule X9(abcc). Therefore, development and medication of such drugs that inhibit TIMP2 outside the cell membrane may be considered to be effective. Similarly, drugs that inhibit MT1-MMP inside the cell membrane might be considered to be useful as well.In addition, upon being clarified the model behaviors we notice the most important pathways in the PWN. They are three parts producingξ91112 (association parts just before the three blue boxes in Figure 5). The PWN may be said to be such a device that continues to produce ξ91112, as explained in the following. In fact, the key molecule X9(abcc) is thus produced, through ξ91112, continuously as far as possible.Figure 5
Looking at PWN with unit variables.The initial concentrationa0 isdispersed in the network, to X4(ab),X7(abc),X9(abcc), X11(abccb), and X12(abccba), according to the first MCL in Table 4. Likewise, b0 and c0 aredispersed in the network, according to latter two MCLs in Table 4, respectively.The units of variables correspond to both the MCLs in Table4 and reactions in Table 5. For example, the first MCL in Table 4 means that the consumption of X1 and production of X4(ab)+X7(abc)+X9(abcc)+X11(abccb)+2X12 are of the same rate in the reaction (a,b) in Table 5.Thus, the initial “mass” continues to exist in the PWN with its components changing to dimer, trimer, tetramer, …: at every moment, such dimer, trimer, …, and hexamercontained in the unit variables coexist. This is true for all the initial masses a0, b0, and c0. By this, in thereactions by unit variables,(1)
X1 versus ξ2581011=X2+X5+X8+2X10+X11,(2)
ξ24=X2+X4 versus ξ3689=X3+2X6+X8+X9,(3)
ξ357=X3+X5+X7 versus ξ357=X3+X5+X7,all possible intermediates of different level oligomers, towards X9, are produced simultaneously at every moment. That is, the PWN system is settled so as to continue to produce X9 as much as possible at every moment.In addition, the node that producesX9(abcc) is not single. There are three such nodes in the PWN, each in reaction of (a,b), (b,c), and (c,c). As in Figure 5, X9 is produced through the production of ξ91112=X9+X11+2X12. Here, this ξ91112 is only produced and no more reused by feedback or additional reactions.Therefore, we may say that the PWN is such a system that producesX9 (through ξ91112) in the triple way (three production nodes) and continues the production as much as possible. This kind of consideration has become possible because the PWN model is analyzed through the unit variables.


### 4.6. On the Simulation Results

We present the time course plots of the theoretical and numerical solutions in Figures6–12. The theoretical solutions above are in good agreement with our simulation results. In Figures 7, 8, 10, and 12, the MCLs (26), (24), (32), and (34), respectively, are confirmed; the solutions are also confirmed. The values of reaction constants and initial values used in the simulations are listed in Table 7. The reaction constants are the same as those used in Watanabe et al. [4], while the initial values are set as follows: the initial values used in [4] are a0=b0=c0=1.0×10-7 [M] based on experimental measurement; if the initial values are adopted as they stand, then the simulation causes a time-delay; that is, the responses become much slower. Therefore, for the sake of convenience in simulation, we take the initial values with the order of 10-6. In addition, in order to show the relationship of the behaviors of solutions and initial values in a more effective way, we set as a0=2.0×10-6 [M]. For example, in Figure 6, we can see that X1(∞)=a0-b0 indeed (see (A.1)). Also, in Figure 8, the MCL can be seen visually.Figure 6
X
1: simulation result and theoretical solution.Figure 7
ξ
2581011: simulation result and theoretical solution and MCL (26).Figure 8
ξ
4791112: simulation result and theoretical solution and MCL (24).Figure 9
ξ
24: simulation result and theoretical solution.Figure 10
ξ
3689: simulation result and theoretical solution and MCL (32).Figure 11
ξ
357: simulation result and theoretical solution.Figure 12
ξ
3689
+
ξ
581011
+
ξ
791112: simulation result and theoretical solution and MCL (36).It should be noted that the simulation results obtained in Watanabe et al. [4] and hence in this paper are based on actual cancer data. The study in this paper provides the exact solutions to the ODE model and hence a mathematically rigorous foundation is added to the simulation results.Concerning the question that if there is a threshold of the key molecule concentration that leads to ECM degradation, we do not know the exact threshold at least presently. It may be considered that, according to the amount of the key molecule, ECM degradation is promoted weakly or strongly. However, as for the unit variableξ3689, it is indicated in Watanabe et al. [4] that, at around b0=100 [nM], transient peak of ξ3689 becomes maximal and the ECM degradation is promoted thereby the most. See Figures S3, S4, and S5 in [4].Figure13 presents the time courses of a simulation result of the key molecule X9 for the model with and without the doubling rule. Here, the initial values of the concentration are a0=2×10-6 [M], b0=2×10-6 [M], and c0=1.8×10-6 [M]. By modifying the model with the doubling rule, the peak at about t=0.25 was increased by approximately 1.5 times.Figure 13
Time course ofX9(t): comparison of the models with or without the doubling rule.Actually, the authors did not begin with the idea of the doubling rules in Section3.2. The only doubling rule we had exploited was the one corresponding to the general basic model in (12). Thus the only coefficients of 2 were on k3X32 and l3X6 in E3, k3X52 and l3X10 in E5, k3X72 and l3X12 in E7, and l2X11 in E11. The solutions to this first version of the model and the modified version in E1–E12 are different, but they have the same initial and asymptotic behaviors. In Figure 13, we show as an example the time series of the key molecule X9(t).Figure14 shows a plot of X9(∞) versus X2(0). As described in Section 3.1, the key molecule of the ECM degradation is abcc (concentration X9), and the amount of abcc that is produced is influenced primarily by b0. If there is too much or too little b0, an insignificant amount of abcc will be produced. As shown in the figure, in both models, X9(∞) has a peak at around b0=5.1×10-7 [M]. The peak increases when the model was modified.Figure 14
X
9
(
∞
) versus X2(0): comparison of the models with or without the doubling rule.


## 4.1. Kinetics of Those Reactions Containinga-bB

According to Table1, the reactions related to a-bB are summarized as (22)X1+X2+X5+X8+2X10+X11⇌l1≡0k1X4+X7+X9+X11+2X12.Looking atX1 and the products (the right hand side) of this reaction, we infer from the flux balance of consumption-production that we may have (23)dX1+X4+X7+X9+X11+2X12dt≡0.By taking summation E1 + E4 + E7 + E9 + E11 + 2 × E12 of the ODEs, we find that (23) is indeed true. Here, X12(abccba) has a coefficient of 2 because the decrease in X1 is equal to the increase in (X4+X7+X9+X11+2X12), since the complex abccba consumes two a molecules, unlike the others: ab,abc,abcc,abccb. Equation (23) is an MCL of form (6a) or (17a). Equation (23) may be considered as an MCL for a0, since (24)X1+X4+X7+X9+X11+2X12t≡const.≡X10=a0by assumption (A2); thus, a local balance holds: a0-X1(t)≡(X4+X7+X9+X11+2X12)(t) for all t≥0.From (22), we also have (25)d-X1+X2+X5+X8+2X10+X11dt≡0,which may be considered to represent an MCL for a0-b0, since (26)-X1+X2+X5+X8+2X10+X11t≡const.≡-X10+X20=-a0+b0.Here, X10 has a coefficient of 2 because the complex bccb has two b- sites for binding with a, which implies a double chance of association, as in reaction (iv). Equation (26) is an MCL of form (6b) or (17b).Note thatE1 and (25) suggest that X1 and (X2+X5+X8+2X10+X11) form a pair of a-bB kinetics, (E1) and (E2581011), as displayed in Table 5. The two ODEs can be solved explicitly using the method for reaction (i) in Section 2; details are shown in Appendix A.1.


## 4.2. Kinetics of Those Reactions ContainingBb-cC

Next, we consider the kinetics of associations/dissociations withBb-cC involved in the ODE system. As can be seen in Table 2, the reactions related to Bb-cC are summarized as (27)X2+X4+X3+2X6+X8+X9⇌l2k2X5+X8+2X10+X11+X7+X9+X11+2X12.In the manner similar to that used for the Aa-bB kinetics, we infer, from the flux balance of consumption-production, that we may have (28)dX2+X5+X8+2X10+X11+X4+X7+X9+X11+2X12dt=0;that is, (29)dX2+X4+X5+X7+X8+X9+2X10+2X11+2X12dt=0.By taking summation E1 + E4 + E7 + E9 + E11 + 2 × E12 of the ODEs, we find that (28) is indeed valid. Equation (28) is an MCL of form (6a) or (17a). Equstion (29) may be viewed as an MCL of b0, since (30)X2+X4+X5+X7+X8+X9+2X10+2X11+2X12t≡const.≡X20=b0.From (27), we also have (31)d-X2-X4+X3+2X6+X8+X9dt=0,which is confirmed by taking the summation −E2 − E4 + E3 + 2 × E6 + E8 + E9 of the ODEs. This is an MCL of form (6b) or (17b). Equation (31) represents an MCL of b0-c0, since (32)-X2-X4+X3+2X6+X8+X9t≡const.≡-X20+X30=-b0+c0.Here, X6 has a coefficient of 2 because it has two c- sites for binding with a or ab.The MCL (31) implies that ξ24(t)≜(X2+X4)(t) and ξ3689(t)≜(X3+2X6+X8+X9)(t) are governed by the same ODE, up to the initial values. From the ODEs E1–E12, it follows that the two variables form a pair of Bb-cC kinetics, (E24) and (E3689), as in Table 4. The variables ξ24(t) and ξ3689(t) can be solved explicitly; details are shown in Appendix A.2.Now, we find that each ofX2(t) and X4(t) can be solved explicitly, as well. This is because the right-hand side of E2 contains only known terms, except for X2: E2′dX2dt=-k1X1t+k2ξ24t+l2-k2b0-c0X2t+l2X1t-a0+b0,by (26) and (32). Hence X2(t) can be calculated by the method of variation of coefficients. See Appendix A.3. In this way, X4(t) can be obtained by X4(t)=ξ24(t)-X2(t).


## 4.3. Kinetics of Those Reactions ContainingCc-cC

We finally consider the kinetics of associations/dissociations ofCc-cC involved in the ODE system. As shown in Table 2, the reactions related to Cc-cC can be summarized as (33)X3+X5+X7+X3+X5+X7⇌l3k32X6+X8+X9+X8+2X10+X11+X9+X11+2X12.In a manner similar to that used previously, we infer from the flux balance of consumption-production that(34)dX3+2X6+X8+X9+X5+X8+2X10+X11+X7+X9+X11+2X12dt=0;that is, (35)dX3+X5+2X6+X7+2X8+2X9+2X10+2X11+2X12dt=0.By taking the summation E3 + E5 + 2 × E6 + E7 + 2 × E8 + 2 × E9 + 2 × E10 + 2 × E11 + 2 × E12 of the ODEs, we find that (35) is indeed valid. Equation (34) is an MCL of form (6a) or (17a). Equation (35) may be viewed as an MCL of c0, since (36)X3+X5+2X6+X7+2X8+2X9+2X10+2X11+2X12t≡const.≡X30=c0.From E3, E5, and E7, we have the ODE (E357), which is displayed in Table 5, and in which we have used (30) and (32). The coefficient of 2 on k3X3(X3+X5+X7)2 in (E357) is a result of the summation of the ODEs E3, E5, and E7. But then, where do the coefficients of 2 on k3(X3+X5+X7) in E3, on k3X5(X3+X5+X7) in E5, and on k3X7(X3+X5+X7) in E7 come from? They stem from the following stoichiometry. For those reactions in Table 3 that correspond to the reaction (ii) in Section 2, (37)X3c+X3c⟷X6cc,X5bc+X5bc⟷X10bccb,X7abc+X7abc⟷X12abccba,the reason is the same as in (12); that is, it is caused by the double rate of consumption/reproduction. Hence, we must count these three reactions with the association/dissociation coefficients on 2k3 and 2l3 in E3, E5, and E7. For the reaction X3(c)+X5(bc)↔X8(bcc) corresponding to reaction (iii), both (38)X3c+X5bc⟷X8bcc,X5bc+X3c⟷X8bccmust be counted. The first reaction in (38) is the one contained in E3, and the second reaction is in E5. In the first reaction, thesubject is X3c, and it associates with theobjectX5(bc). We will name this molecule c=c1. In the second reaction, thesubject is X5bc, and it associates with theobjectX3(c). We will name this molecule c=c2. In other words, a molecule of c (=object) is associated with a molecule of bc (=subject). In general, c=c1 and c2 are different molecules.The first reaction must be counted as a reaction ofX3c. Also, seen as a reaction of X5bc, the second reaction must also be counted. Both reactions take place at the same time. The association/dissociation coefficients are k3 and l3, respectively, for both reactions.However, now consider, for example, howX3(c) is consumed doubly by the amount of c1 and c2. In Table 3, other combinations of this kind are (39)X5bc+X7abc⟷X11abccb,X7abc+X5bc⟷X11abccb,X3c+X7abc⟷X9abcc,X7abc+X3c⟷X9abcc.The contribution of the reactions in (38) and (39) are therefore doubled. Combining the doubling of the reactions in (38), we obtain the third reaction in Table 5, which causes the term -2k3(X3+X5+X7)2 in E357 and the ODEs E3, E5, and E7.Finally, in consideringX3(c), the dissociation term l3(X8+X9) need not be doubled, since only one molecule of c is produced from a molecule of bcc: bc-c; also, this dissociation is not involved in E5.By the above arguments for the doubling rules, we have obtained ingeniously well-formed MCLs and reaction laws, as shown in Table4, respectively. This may suggest that our doubling rules are justified.Also, the above argument ofCc-cC kinetics may imply the following: suppose that, for association of molecules, p1,p2,…,pn, say, some combinations of coupling pi-pj (1≤i,j≤n), are possible and others not. Then, it does not hold that (40)dX1+⋯+Xndt=-2kX1+⋯+Xn2+dissociationtermin general. However, if all of pi has a common binding cite, that is, all pi are of the form -cC, then we do have (40). In the latter case, the doubling rule holds for (X1+⋯+Xn)as a unit, as if (p1,p2,…,pn) are a single molecule with binding cite -c, and the elementary doubling rule (12) holds with [A] in (12) replaced by (X1+⋯+Xn).


## 4.4. Solutions to Remaining Variables

Solutions to the remaining variables can also be obtained.X3 and X5 can be found by the method of variation of parameters, as in the case of X2. X7 can be found by X7(t)=ξ357(t)-(X3+X5)(t), where ξ357(t) is the solution to (X3+X5+X7)(t). By applying the solution ξ357 to E3, we obtain the solution to X3(t) as E3′dX3dt=-k2-l3ξ24+2k3-l2ξ357+l2+l3X3t+l2ξ357+l3ξ24-l3b0-c0by (30) and (32); note that ξ24 and ξ357 have been obtained already. Similarly, X6, X8, and X10 can also be found by the method of variation of coefficient.All ofX1 to X8 and X10 have thus been obtained in principle. The remaining ones, X9, X11, and X12, can be obtained by solving (41)X9t=ξ3689t-X3+2X6+X8t,X11t=X1t-X2+X5+X8+2X10t-a0+b0,X12t=12a0-X1+X4+X7+X9+X11t(see (A.8), (26), and (24), resp.). The processes for obtaining these solutions are listed in Table 6.


## 4.5. Biomedical Implication

As described in Section3.1, an appropriate amount of X2(b) (TIMP2) is necessary to obtain a substantial amount of the key molecule X9(abcc). Therefore, development and medication of such drugs that inhibit TIMP2 outside the cell membrane may be considered to be effective. Similarly, drugs that inhibit MT1-MMP inside the cell membrane might be considered to be useful as well.In addition, upon being clarified the model behaviors we notice the most important pathways in the PWN. They are three parts producingξ91112 (association parts just before the three blue boxes in Figure 5). The PWN may be said to be such a device that continues to produce ξ91112, as explained in the following. In fact, the key molecule X9(abcc) is thus produced, through ξ91112, continuously as far as possible.Figure 5
Looking at PWN with unit variables.The initial concentrationa0 isdispersed in the network, to X4(ab),X7(abc),X9(abcc), X11(abccb), and X12(abccba), according to the first MCL in Table 4. Likewise, b0 and c0 aredispersed in the network, according to latter two MCLs in Table 4, respectively.The units of variables correspond to both the MCLs in Table4 and reactions in Table 5. For example, the first MCL in Table 4 means that the consumption of X1 and production of X4(ab)+X7(abc)+X9(abcc)+X11(abccb)+2X12 are of the same rate in the reaction (a,b) in Table 5.Thus, the initial “mass” continues to exist in the PWN with its components changing to dimer, trimer, tetramer, …: at every moment, such dimer, trimer, …, and hexamercontained in the unit variables coexist. This is true for all the initial masses a0, b0, and c0. By this, in thereactions by unit variables,(1)
X1 versus ξ2581011=X2+X5+X8+2X10+X11,(2)
ξ24=X2+X4 versus ξ3689=X3+2X6+X8+X9,(3)
ξ357=X3+X5+X7 versus ξ357=X3+X5+X7,all possible intermediates of different level oligomers, towards X9, are produced simultaneously at every moment. That is, the PWN system is settled so as to continue to produce X9 as much as possible at every moment.In addition, the node that producesX9(abcc) is not single. There are three such nodes in the PWN, each in reaction of (a,b), (b,c), and (c,c). As in Figure 5, X9 is produced through the production of ξ91112=X9+X11+2X12. Here, this ξ91112 is only produced and no more reused by feedback or additional reactions.Therefore, we may say that the PWN is such a system that producesX9 (through ξ91112) in the triple way (three production nodes) and continues the production as much as possible. This kind of consideration has become possible because the PWN model is analyzed through the unit variables.


## 4.6. On the Simulation Results

We present the time course plots of the theoretical and numerical solutions in Figures6–12. The theoretical solutions above are in good agreement with our simulation results. In Figures 7, 8, 10, and 12, the MCLs (26), (24), (32), and (34), respectively, are confirmed; the solutions are also confirmed. The values of reaction constants and initial values used in the simulations are listed in Table 7. The reaction constants are the same as those used in Watanabe et al. [4], while the initial values are set as follows: the initial values used in [4] are a0=b0=c0=1.0×10-7 [M] based on experimental measurement; if the initial values are adopted as they stand, then the simulation causes a time-delay; that is, the responses become much slower. Therefore, for the sake of convenience in simulation, we take the initial values with the order of 10-6. In addition, in order to show the relationship of the behaviors of solutions and initial values in a more effective way, we set as a0=2.0×10-6 [M]. For example, in Figure 6, we can see that X1(∞)=a0-b0 indeed (see (A.1)). Also, in Figure 8, the MCL can be seen visually.Figure 6
X
1: simulation result and theoretical solution.Figure 7
ξ
2581011: simulation result and theoretical solution and MCL (26).Figure 8
ξ
4791112: simulation result and theoretical solution and MCL (24).Figure 9
ξ
24: simulation result and theoretical solution.Figure 10
ξ
3689: simulation result and theoretical solution and MCL (32).Figure 11
ξ
357: simulation result and theoretical solution.Figure 12
ξ
3689
+
ξ
581011
+
ξ
791112: simulation result and theoretical solution and MCL (36).It should be noted that the simulation results obtained in Watanabe et al. [4] and hence in this paper are based on actual cancer data. The study in this paper provides the exact solutions to the ODE model and hence a mathematically rigorous foundation is added to the simulation results.Concerning the question that if there is a threshold of the key molecule concentration that leads to ECM degradation, we do not know the exact threshold at least presently. It may be considered that, according to the amount of the key molecule, ECM degradation is promoted weakly or strongly. However, as for the unit variableξ3689, it is indicated in Watanabe et al. [4] that, at around b0=100 [nM], transient peak of ξ3689 becomes maximal and the ECM degradation is promoted thereby the most. See Figures S3, S4, and S5 in [4].Figure13 presents the time courses of a simulation result of the key molecule X9 for the model with and without the doubling rule. Here, the initial values of the concentration are a0=2×10-6 [M], b0=2×10-6 [M], and c0=1.8×10-6 [M]. By modifying the model with the doubling rule, the peak at about t=0.25 was increased by approximately 1.5 times.Figure 13
Time course ofX9(t): comparison of the models with or without the doubling rule.Actually, the authors did not begin with the idea of the doubling rules in Section3.2. The only doubling rule we had exploited was the one corresponding to the general basic model in (12). Thus the only coefficients of 2 were on k3X32 and l3X6 in E3, k3X52 and l3X10 in E5, k3X72 and l3X12 in E7, and l2X11 in E11. The solutions to this first version of the model and the modified version in E1–E12 are different, but they have the same initial and asymptotic behaviors. In Figure 13, we show as an example the time series of the key molecule X9(t).Figure14 shows a plot of X9(∞) versus X2(0). As described in Section 3.1, the key molecule of the ECM degradation is abcc (concentration X9), and the amount of abcc that is produced is influenced primarily by b0. If there is too much or too little b0, an insignificant amount of abcc will be produced. As shown in the figure, in both models, X9(∞) has a peak at around b0=5.1×10-7 [M]. The peak increases when the model was modified.Figure 14
X
9
(
∞
) versus X2(0): comparison of the models with or without the doubling rule.


## 5. Concluding Remarks

We presented a method for solving a nonlinear ODE system from cell biology. With the setting in this paper, the ODEs become a completely integrable system so that they can be solved explicitly. The key idea was to use the MCL to obtain linear relations of the ODE variables that would be valid for allt≥0. The previous analysis of such ODE systems was primarily based on the balance of flux at equilibrium (t=∞) or by computer simulation (for t<∞). The theoretical solutions to the ODEs are in complete agreement with those obtained by simulation.We gained several new insights in analyzing and modeling the PWNs: in the MALs, a coefficient of 2 must be attached sometimes to appropriate molecules. Determining all of the doubling laws is crucial for obtaining an integrable system of ODEs. We have elucidated some situations in which the doubling law is required. We also note that the total PWN can be grouped into several local PWNs in which local linear relationships hold. This may suggest a way to obtain an appropriate view of the composition of network motifs in systems biology.In a future study, we will try to determine a sufficient condition such that the given ODE system is completely integrable. It may be desirable to clarify under which conditions the ODE system framework is sufficiently “close” to a realistic PDE system; this would be useful for the development of therapies or drugs.



---

*Source: 5924270-2017-05-30.xml*