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

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 called doubling 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.


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 all ≥ 0. As a result, the nonlinear ODEs presented by the new units of variables turns out to be a completely 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 all ≥ 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 solved as 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 including doubling rule is necessary as in the following, which has not been considered properly so far.  Dissoc. const.
The organization of the paper is as follows. Section 2 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, the doubling 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.

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 wellknown 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 by and , respectively.
(i) Association of two different molecules/dissociation of their product: (ii) Association of the same molecules/dissociation of their product: (iii) Association with a modified molecule/dissociation of their product: (iv) Association with a homodimer/dissociation of their product: Note that + is the same as reaction (iii).
( , ): (X 3 + X 5 + X 7 ) + ( 3 + 5 + 7 ) 3 3 (2X 6 + X 8 + X 9 + X 8 + 2X 10 + X 11 + X 9 + 11 + 2X 12 ) ⇐⇒ where the last expression is the factorization of the quadratic polynomial in the second line. If 1 ̸ = 2 , this leads to , and hence the explicit solution for ≥ 0. Here, 1 , 2 are In the first equation of (12), the coefficient of 2 stems from the fact that, in the reaction + → , two molecules of are consumed; similarly, in the dissociation, two molecules of are produced. This is called a double rate of consumption/reproduction. The ODEs (12) imply the MCL for ≥ 0. Here 1 , 2 are The ODEs (16) can be solved in a manner similar to that used for reactions (i) and (ii). In fact, the association ofreduces to reaction (i), and the association of -reduces to reaction (ii).
Finally, for reaction (iv), we have the following MALs and the MCL: The coefficient of 2 in (17) arises because a molecule of may bind to either of the two molecules. The solution to these ODEs (17) is the same as that for reaction (i) but with replaced by 2 . This is called a double 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 equations as groups to those of the elementary reaction processes.

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.

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 by , , and , 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) Although (MT1-MMP) is an ECM degradation enzyme as well, it is usually inactivated by (TIMP2) immediately after vesicle transportation to the cell membranes. However, activates (MMP2), which is another ECM degradation enzyme, as follows.
(1) The molecules that penetrate the cell membrane form homodimers on the membrane.
(2) On one of the in the dimer --, the heterodimer (pro-MMP2, TIMP2) is coupled to produce the tetramer ; alternatively, -may be coupled with the heterotrimer -to produce the tetramer. Thus, is associated with via .
(3) One of the in that is not coupled with cuts the connection between and .
(4) that was cut and released then becomes the activated one.
This process is thus the associations/dissociations of the three molecules , , and , and is the key molecule in the activation of . In order to produce certain amount of activated (MMP2), there must be sufficient . However, if there is too much , then the remaining vacant site of tends to associate with -, which results in the fact that or will be produced preferentially to , and an insignificant amount of will be activated. See Hoshino et al. [3] and then Watanabe et al. [4]. The molecular binding rules that follow from this scenario are that and do not couple with each other, but with , and may form the homodimer (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) has a site of binding only with .
(2) has a site of binding with and another site of binding with . , , and . The PWN of association or dissociation of these complexes are shown in Figure 3.
Thus, a quantitative model of ECM degradation can be based on the kinetics of associations/dissociations of the three proteins ( : MMP2; : TIMP2; and : 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 , , and 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 ( , , and ) 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; see (E1)-(E12). Now, according to the binding rules discussed above, we have the following three associations or dissociations in the PWN: (I) Association of -(association constant 1 ): The molecules and do not dissociate after they are once associated. Thus we set 1 ≡ 0.
(II) Association/dissociation of -(association constant 2 ; dissociation constant 2 ): (III) Association/dissociation of -(association constant 3 ; dissociation constant 3 ): For reactions (I)-(III), there are twelve molecules in the PWN in Figure 3: , , , , , , , . We denote their concentrations as 1 = [ ], 2 = [ ], 3 = [ ], and so on. The three association/dissociation groups are presented in Tables 1, 2, and 3. We make the following basic assumptions: ( 1) For the association/dissociation of the modified molecules and , we use the same association/dissociation constants as those used for the elementary reaction -. The same is assumed forand -. In order to obtain the correct ODEs, it is necessary to extract the appropriate doubling rules for the reaction coefficients in the PWN. We will explain it in the next subsection. Figure 3 and Tables 1-3, we can write down the following twelve second-order nonlinear ODEs:
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.
Here, X 10 has a coefficient of 2 because the complex bccb has two bsites for binding with , which implies a double chance of association, as in reaction (iv). Equation (26) is an MCL of form (6b) or (17b). Note that (E1) and (25) suggest that 1 and (X 2 +X 5 +X 8 + 2X 10 + X 11 ) form a pair of -kinetics, ( 1 ) and ( 2581011 ), 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.
In other words, a molecule of (=object) is associated with a molecule of bc (=subject). In general, c = 1 and 2 are different molecules. The first reaction must be counted as a reaction of X 3 (c). Also, seen as a reaction of X 5 (bc), the second reaction must also be counted. Both reactions take place at the same time. The association/dissociation coefficients are 3 and 3 , respectively, for both reactions.
Finally, in considering 3 ( ), the dissociation term 3 ( 8 + 9 ) need not be doubled, since only one molecule of is produced from a molecule of : -; also, this dissociation is not involved in (E5).
4.5. Biomedical Implication. As described in Section 3.1, an appropriate amount of 2 ( ) (TIMP2) is necessary to obtain a substantial amount of the key molecule 9 ( ). 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 9 ( ) is thus produced, through 91112 , continuously as far as possible.
The initial concentration 0 is dispersed in the network, to X 4 (ab), X 7 (abc), X 9 (abcc), X 11 (abccb), and X 12 (abccba), according to the first MCL in Table 4. Likewise, 0 and 0 are dispersed in the network, according to latter two MCLs in Table 4, respectively.
The units of variables correspond to both the MCLs in Table 4 and reactions in Table 5. For example, the first MCL in Table 4 means that the consumption of 1 and production of X 4 (ab) + X 7 (abc) + X 9 (abcc) + X 11 (abccb) + 2X 12 are of the same rate in the reaction ( , ) in Table 5.
Therefore, we may say that the PWN is such a system that produces 9 (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 Figures  6-12. The theoretical solutions above are in good agreement with our simulation results. In Figures 7, 8, 10, and 12, the Sum of the two X 1 (t) (X 4 + X 7 + X 9 + X 11 + 2X 12 )(t)  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 0 = 0 = 0 = 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 0 = 2.0 × 10 −6 [M]. For example, in Figure 6, we can see that

× 10 −6
Sum of the two b 0 − (X 2 + X 4 )(t) (X 3 + X 6 + 2X 8 + X 9 )(t) X 3 + X 5 + X 7 (t), numerical X 3 + X 5 + X 7 (t), theoretical 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 0 = 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].  6 8 1 0 2 t (X 3 + X 6 + 2X 8 + X 9 )(t) (X 7 + X 9 + X 11 + 2X 12 )(t) Sum of the three (X 5 + X 8 + 2X 10 + X 11 )(t) By modifying the model with the doubling rule, the peak at about = 0.25 was increased by approximately 1.5 times. Actually, the authors did not begin with the idea of the doubling rules in Section 3.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 3 2 3 and 3 6 in (E3), 3 2 5 and 3 10 in (E5), 3 2 7 and 3 12 in (E7), and 2 11 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 9 ( ). Figure 14 shows a plot of 9 (∞) versus 2 (0). As described in Section 3.1, the key molecule of the ECM degradation is (concentration 9 ), and the amount of that is produced is influenced primarily by 0 . If there is too much or too little 0 , an insignificant amount of will be produced. As shown in the figure, in both models, 9 (∞) has a peak at around 0 = 5.1 × 10 −7 [M]. The peak increases when the model was modified.

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 all ≥ 0. The previous analysis of such ODE systems was primarily based on the balance of flux at equilibrium ( = ∞) or by computer simulation (for < ∞). 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.