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 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.
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 all . 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 . 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.
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.
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 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).
We now show that the ODEs for the MALs for reactions (i)–(iv) can be solved explicitly. First, the MALs for reaction (i) areThe resulting system of ODEs can be solved explicitly: (5) implies the MCLsfor some positive constants and , which are typically the initial values of and , respectively. Then, substituting into the third equation in (5), we obtain an equation in only the variable : where the last expression is the factorization of the quadratic polynomial in the second line. If , this leads to or , and hence the explicit solution for . Here, are If and are distinct and so that in (10), then we find that is monotonically decreasing and nonnegative. Solutions to and 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 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 a positive constant . Then, substituting into the first equation of (12), we have for . Here are Since , we have . As in case (i), the solution is monotonically decreasing and nonnegative. The solution to can also be obtained by the MCL (13).
For reaction (iii), we have the MALs and MCL as shown in the following:The ODEs (16) can be solved in a manner similar to that used for reactions (i) and (ii). In fact, the association of - reduces 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 arises because a molecule of may bind to either of the two molecules. The solution to these ODEs is the same as that for reaction (i) but with replaced by . 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.
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 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 .(3) has a site of binding with and another site of binding with itself.(4)There are no other binding patterns.
By inspection, we see that there appear to be nine complexes (dimers to hexamer) obtainable from the monomers , , 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 –.
Now, according to the binding rules discussed above, we have the following three associations or dissociations in the PWN:(I)Association of - (association constant ): The molecules and do not dissociate after they are once associated. Thus we set .(II)Association/dissociation of - (association constant ; dissociation constant ): (III)Association/dissociation of - (association constant ; dissociation constant ):
For reactions (I)–(III), there are twelve molecules in the PWN in Figure 3: , , , , , , , , , , , and . We denote their concentrations as , , , and so on. The three association/dissociation groups are presented in Tables 1, 2, and 3. We make the following basic assumptions: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 for - and -.The initial concentrations for are all : , .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.
3.2. Pathway Network Dynamics and Doubling Rules
From Figure 3 and Tables 1–3, we can write down the following twelve second-order nonlinear ODEs:
Note that ODEs – incorporate doubling rules. The coefficient of 2 on(1) in ,(2) in ,(3) in corresponds to the double rate of consumption/reproduction in reaction (ii). See Figure 4. The coefficient of 2 on (1) in ,(2) in ,(3) in ,(4) in ,(5) in ,(6) in ,which corresponds to - 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.
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 Containing -
According to Table 1, the reactions related to - are summarized as
Looking at and the products (the right hand side) of this reaction, we infer from the flux balance of consumption-production that we may have By taking summation + + + + + 2 × of the ODEs, we find that (23) is indeed true. Here, has a coefficient of 2 because the decrease in is equal to the increase in , since the complex consumes two molecules, unlike the others: . Equation (23) is an MCL of form (6a) or (17a). Equation (23) may be considered as an MCL for , since by assumption (2); thus, a local balance holds: for all .
From (22), we also have which may be considered to represent an MCL for , since Here, has a coefficient of 2 because the complex has two - sites 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 and (25) suggest that and form a pair of - kinetics, and , 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 Containing -
Next, we consider the kinetics of associations/dissociations with - involved in the ODE system. As can be seen in Table 2, the reactions related to - are summarized as In the manner similar to that used for the - kinetics, we infer, from the flux balance of consumption-production, that we may have that is, By taking summation + + + + + 2 × 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 , since
From (27), we also have which is confirmed by taking the summation − − + + × + + of the ODEs. This is an MCL of form (6b) or (17b). Equation (31) represents an MCL of , since Here, has a coefficient of 2 because it has two - sites for binding with or .
The MCL (31) implies that and are governed by the same ODE, up to the initial values. From the ODEs –, it follows that the two variables form a pair of - kinetics, and , as in Table 4. The variables and can be solved explicitly; details are shown in Appendix A.2.
Now, we find that each of and can be solved explicitly, as well. This is because the right-hand side of contains only known terms, except for : by (26) and (32). Hence can be calculated by the method of variation of coefficients. See Appendix A.3. In this way, can be obtained by .
4.3. Kinetics of Those Reactions Containing -
We finally consider the kinetics of associations/dissociations of - involved in the ODE system. As shown in Table 2, the reactions related to - can be summarized as
In a manner similar to that used previously, we infer from the flux balance of consumption-production that that is, By taking the summation + + 2 × + + 2 × + 2 × + 2 × + 2 × + 2 × 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 , since From , , and , we have the ODE , which is displayed in Table 5, and in which we have used (30) and (32). The coefficient of 2 on in is a result of the summation of the ODEs , , and . But then, where do the coefficients of 2 on in , on in , and on in come from? They stem from the following stoichiometry. For those reactions in Table 3 that correspond to the reaction (ii) in Section 2, 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 and in , , and . For the reaction corresponding to reaction (iii), both must be counted. The first reaction in (38) is the one contained in , and the second reaction is in . In the first reaction, the subject is , and it associates with the object . We will name this molecule . In the second reaction, the subject is , and it associates with the object . We will name this molecule . In other words, a molecule of (=object) is associated with a molecule of (=subject). In general, and are different molecules.
The first reaction must be counted as a reaction of . Also, seen as a reaction of , the second reaction must also be counted. Both reactions take place at the same time. The association/dissociation coefficients are and , respectively, for both reactions.
However, now consider, for example, how is consumed doubly by the amount of and . In Table 3, other combinations of this kind are 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 in and the ODEs , , and .
Finally, in considering , the dissociation term need not be doubled, since only one molecule of is produced from a molecule of : -; also, this dissociation is not involved in .
By the above arguments for the doubling rules, we have obtained ingeniously well-formed MCLs and reaction laws, as shown in Table 4, respectively. This may suggest that our doubling rules are justified.
Also, the above argument of - kinetics may imply the following: suppose that, for association of molecules, , say, some combinations of coupling - (), are possible and others not. Then, it does not hold that in general. However, if all of has a common binding cite, that is, all are of the form -, then we do have (40). In the latter case, the doubling rule holds for as a unit, as if are a single molecule with binding cite -, and the elementary doubling rule (12) holds with in (12) replaced by .
4.4. Solutions to Remaining Variables
Solutions to the remaining variables can also be obtained. and can be found by the method of variation of parameters, as in the case of . can be found by , where is the solution to . By applying the solution to , we obtain the solution to as by (30) and (32); note that and have been obtained already. Similarly, , , and can also be found by the method of variation of coefficient.
All of to and have thus been obtained in principle. The remaining ones, , , and , can be obtained by solving (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 Section 3.1, an appropriate amount of (TIMP2) is necessary to obtain a substantial amount of the key molecule . 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 (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 , as explained in the following. In fact, the key molecule is thus produced, through , continuously as far as possible.
The initial concentration is dispersed in the network, to , , and , according to the first MCL in Table 4. Likewise, and 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 and production of are of the same rate in the reaction 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 hexamer contained in the unit variables coexist. This is true for all the initial masses , , and . By this, in the reactions by unit variables,(1) versus ,(2) versus ,(3) versus ,all possible intermediates of different level oligomers, towards , are produced simultaneously at every moment. That is, the PWN system is settled so as to continue to produce as much as possible at every moment.
In addition, the node that produces is not single. There are three such nodes in the PWN, each in reaction of , , and . As in Figure 5, is produced through the production of . Here, this is only produced and no more reused by feedback or additional reactions.
Therefore, we may say that the PWN is such a system that produces (through ) 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 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 [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 . In addition, in order to show the relationship of the behaviors of solutions and initial values in a more effective way, we set as [M]. For example, in Figure 6, we can see that indeed (see (A.1)). Also, in Figure 8, the MCL can be seen visually.
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 , it is indicated in Watanabe et al. [4] that, at around [nM], transient peak of becomes maximal and the ECM degradation is promoted thereby the most. See Figures S3, S4, and S5 in [4].
Figure 13 presents the time courses of a simulation result of the key molecule for the model with and without the doubling rule. Here, the initial values of the concentration are [M], [M], and [M]. By modifying the model with the doubling rule, the peak at about 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 and in , and in , and in , and in . The solutions to this first version of the model and the modified version in – 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 .
Figure 14 shows a plot of versus . As described in Section 3.1, the key molecule of the ECM degradation is (concentration ), and the amount of that is produced is influenced primarily by . If there is too much or too little , an insignificant amount of will be produced. As shown in the figure, in both models, has a peak at around [M]. The peak increases when the model was modified.
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 all . 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.
Appendix
A. Explicit Solutions to the Unit Variables
The solutions to – and related materials are presented as follows.
A.1. Solution to : Single Riccati Equation
Theorem A.1. One has the following solution to :
Corollary A.2. is decreasing monotonically:
Thus is nonnegative as well. Let us denote .
Corollary A.3. One has the following solution to :
Thus is decreasing monotonically: and is nonnegative.
A.2. Solution to : Single Riccati Equation
Theorem A.4. One has the following solution to : where
We note that . Since , we have the following.
Corollary A.5. is decreasing monotonically: Thus is nonnegative as well.
Corollary A.6. One has the following solution to : Thus is decreasing monotonically: and is nonnegative.
The solutions and are thus obtained irrespective of or . Each of the solutions to and , however, does depend on , or as seen in the following.
A.3. Solution to and Thus : Method of Variation of Constant
Let , , and .
Theorem A.7. One hasrespectively. Here is given by
The solutions turn out to be nonnegative in either of the three cases. Now is given by , although we do not write down the solution here because it is complicated like .
A.4. Solution to : Single Riccati Equation
The solution to is obtained, similarly as . whererespectively. Properties of positivity and monotone decreasing of are valid as well.
The solutions to other variables, , , , , , , , , and , are too complicated to write down and omitted.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The first and the last authors are supported by JSPS KAKENHI Grant-in-Aid for Exploratory Research, no. 24654023. The last author is supported also by JST-CREST and KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas, no. 16H06576. The authors would like to express their appreciation.