Control and inhibition analysis of complex formation processes

Abstract

Background

Proteolytic degradation of the extracellular matrix (ECM) is a key event in tumour metastasis and invasion. Matrix metalloproteinases (MMPs) are a family of endopeptidases that degrade most of the components of the ECM. Several broad-spectrum MMP inhibitors (MMPIs) have been developed, but have had little success due to side effects. Thus, it is important to develop mathematical methods to provide new drug treatment strategies. Matrix metalloproteinase 2 (MMP2) activation occurs via a mechanism involving complex formation that consists of membrane type 1 MMP (MT1-MMP), tissue inhibitor of matrix metalloproteinase 2 (TIMP2) and MMP2. Here, we focus on developing a method for analysing the complex formation process.

Results

We used control analysis to investigate inhibitor responses in complex formation processes. The essence of the analysis is to define the response coefficient which measures the inhibitory efficiency, a small fractional change of concentration of a targeting molecule in response to a small fractional change of concentration of an inhibitor. First, by using the response coefficient, we investigated models for general classes of complex formation processes: chain reaction systems composed of ordered steps, and chain reaction systems and site-binding reaction systems composed of unordered multi-branched steps. By analysing the ordered step models, we showed that parameter-independent inequalities between the response coefficients held. For the unordered multi-branched step models, we showed that independence of the response coefficients with respect to equilibrium constants held. As an application of our analysis, we discuss a mathematical model for the MMP2 activation process. By putting the experimentally derived parameter values into the model, we were able to conclude that the TIMP2 and MMP2 interaction is the most efficient interaction to consider in selecting inhibitors.

Conclusions

Our result identifies a new drug target in the process of the MMP2 activation. Thus, our analysis will provide new insight into the design of more efficient drug strategies for cancer treatment.

Background

Metastasis and invasion are a major cause of death in cancer patients, and thus preventing this secondary spread of the tumour is an important aspect of cancer therapy. A prerequisite for the migration of endothelial cells through the extracellular matrix (ECM) is the initiation of a biochemical pathway responsible for the proteolytic degradation of these structural barriers. Matrix metalloproteinases (MMPs) are a family of endopeptidases that degrade most of the components of the ECM [1, 2]. Several broad-spectrum MMP inhibitors (MMPIs) have been developed, some of which have been used in clinical trials for cancer treatment [3]. However, these MMP inhibitors have had little success due to side effects, implying that these clearly lacked selectivity in their action. Most MMPs are inhibited by MMPIs that bind to the active sites of MMPs. Similarities in active sites of MMPs pose obstacles to the design of specific inhibitors. It is thus important to find alternatives to these approaches to increase specificity. Therefore, it is worth developing mathematical methods in analysing biochemical reaction pathways to provide new drug treatment strategies.

Matrix metalloproteinase 2 (MMP2) was proposed as a potential therapeutic target, based on its high-level expression in many human tumours and its ability to degrade type IV collagen [4]. The activation of MMP2 proenzyme is processed by membrane type 1 matrix metalloproteinase (MT1-MMP) [5–7]. Under physiologic conditions, MMP2 is secreted as a latent form, pro-MMP2, and it has been established that its activation occurs via a mechanism involving a complex formation that consists of MT1-MMP, tissue inhibitor of matrix metalloproteinase 2 (TIMP2) and pro-MMP2. Here, focusing on the MMP2 activation process, we develop a mathematical method which quantifies a response of systems to an inhibitor and classifies interactions in order of the inhibitory efficiency. The goal of this study is to identify the most efficient interaction to consider in selecting inhibitors.

In order to achieve this goal, we use control analysis, a method of sensitivity analysis which is widely used in many fields (see for example [8]). Methods of obtaining quantitative measures of control in metabolic pathways were developed by Kacser and Burns [9] and Heinrich and Rapoport [10, 11] (for a review, see [12]). The essence of the analysis is to define the response coefficient which measures the inhibitory efficiency, a small fractional change of concentration of a targeting molecule in response to a small fractional change of concentration of an inhibitor. First, using the response coefficient, we investigate models for general classes of complex formation processes: chain reaction systems composed of ordered steps, and chain reaction systems and site-binding reaction systems composed of unordered multi-branched steps. By analysing the ordered step models, we show that parameter-independent inequalities between response coefficients hold. In the unordered multi-branched step models, we assume there are no cooperative reactions, i.e. the equilibrium constant for binding a site of molecule B to a site of another molecule A is independent of whether any of the other sites of molecule A are occupied. This assumption satisfies the detailed balance condition, but it puts a stronger constraint on the system (the importance of the detailed balance condition in kinetic modelling was discussed in [13–15]). Under this assumption, we show that independence of the response coefficients with respect to the equilibrium constants holds. These results indicate that the inhibitory efficiency depends on the topology of pathway networks. As an application of our analysis, we investigate a mathematical model for the MMP2 activation process [16]. In the model, MT1-MMP is dimerized, bound to TIMP2 and forms a quadruple complex by binding proMMP2, MT1-MT1T2M2, which contributes to activation of proMMP2. We try to identify the most efficient interaction to consider in selecting inhibitors by putting the experimentally derived parameter values into the model.

Therefore, our method can serve as a tool for quantifying the inhibitory efficiency and allows us to determine the most efficient method of selecting the inhibitor in complex formation processes. Application of the method to the MMP model may provide a new tool for designing more efficient drug strategies for cancer treatment.

Results

The response coefficient as a measure of the inhibitory efficiency

Here, we explain fundamentals of control analysis. The key point of the analysis is to introduce the response coefficient, which describes how a variable such as molecular concentration responds to variation of parameters such as rate constants. In this paper, we would like to consider the response of an inhibitor concentration. Mathematically, the response is expressed as

$$\begin{array}{l}R=\underset{\delta \left[I\right]\to 0}{lim}\frac{\delta F/F}{\delta \left[I\right]/\left[I\right]}\\ =\frac{\partial F}{\partial \left[I\right]}\frac{\left[I\right]}{F}=\frac{\partial lnF}{\partial ln\left[I\right]}\text{,}\end{array}$$

where F denotes the objective function, which will be taken as a molecular concentration in the steady state of the system. Therefore, the response coefficient R represents a fractional change of the objective function F, δF/F, in response to a fractional change in an inhibitor concentration, δ[I]/[I], in the limit as δ[I] tends to zero. Introducing an inhibitor requires construction of a new model which consists of the original species and a new inhibitor. Instead of introducing new inhibitors, we here consider the change of the objective function in response to changes in the equilibrium constants. To see why this alternative method is valid, let us consider a simple Michaelis–Menten type enzymatic reaction

$$\text{E}+\text{S}\begin{array}{c}\hfill \to \hfill \\ \hfill \leftarrow \hfill \end{array}\text{ES}\to E+P\text{.}$$

Suppression of concentration [ES] leads to suppression of P production. Hence, to make analysis simpler, we consider a simpler reaction scheme

$$\text{E}+\text{S}\begin{array}{c}\hfill \to \hfill \\ \hfill \leftarrow \hfill \end{array}\text{ES}\text{.}$$

Disturbance of this reaction by an inhibitor in a competitive way is expressed as

$$\begin{array}{l}\text{E}+\text{S}\begin{array}{c}\hfill \to \hfill \\ \hfill \leftarrow \hfill \end{array}\text{ES}\\ +\\ \text{I ,}\\ \uparrow \downarrow \\ \text{EI}\end{array}$$

The steady state solution of this reaction system can be written as

$$\begin{array}{l}\left[ES\right]=\frac{\left[S\right]}{\alpha K+\left[S\right]}{\left[E\right]}_T,\\ \left[E\right]+\left[EI\right]=\frac{\alpha K}{\alpha K+\left[S\right]}{\left[E\right]}_{T\text{,}}\end{array}$$

where [E] _T is the total concentration of the enzyme. The concentration [ES] is the enzyme concentration which contributes to P production, while the concentration [E] + [EI] is the enzyme concentration which does not contribute to P production. The parameter K, defined by K = kd/ka, is the equilibrium constant of the enzyme and substrate binding reaction, where kd is the dissociation rate constant and ka is the association rate constant. The coefficient α is written as $\alpha =1+\frac{\left[I\right]}{K_I}$, where K _I is the equilibrium constant of the enzyme and inhibitor binding reaction. Because α is always greater than 1, an addition of the inhibitor is equivalent to an increase of the equilibrium constant K. Therefore, in this paper, instead of adding inhibitors, we perform an inhibition analysis by changing equilibrium constants.

Let us take a look at the basic properties of the response coefficient for a simple model. The response coefficient for F with respect to K is written as

$$R=\frac{\partial F}{F}/\frac{\partial K}{K}=\frac{\partial lnF}{\partial lnK}\text{.}$$

(1)

At the steady state, F is a general function of total concentrations of molecules and equilibrium constants, and thus we can write F = F([E] _T , [S] _T , K). Let us see how the response coefficient behaves as the value of the equilibrium constant changes. First, suppose that the total concentration of the substrate is much greater than that of the enzyme, i.e. [E] _T  < < [S] _T . In this case, the system can be reduced to a linear system. The steady state solution is written simply as

$$\left[ES\right]=\frac{{\left[S\right]}_T}{K+{\left[S\right]}_T}{\left[E\right]}_T\text{.}$$

The response coefficient for [ES] with respect to K is calculated as

$$R^{\mathit{ES}}=\frac{\partial ln\left[ES\right]}{\partial lnK}=-\frac{K}{K+{\left[S\right]}_T}\text{.}$$

The minus sign means that a decrease of the complex concentration [ES] is caused by an increase of K. The response coefficient R^ES goes to 0 as K goes to 0, while R^ES goes to −1 as K goes to ∞. This indicates that the inhibition becomes more efficient as K gets larger. Next, to explore behaviours of the response coefficient in the nonlinear regime (the enzyme concentration is in the same order as the substrate concentration) of the system, we use numerical simulations. We calculate R^ES for three different total substrate concentrations, [S] _T  = 1 nM, 100 nM and 10000 nM (Figure 1A). All three curves behave in the same manner in that R^ES goes to 0 as K goes to 0, while R^ES goes to −1 as K becomes large, as described above. To compare these three curves quantitatively, we introduce the Hill coefficient n _H , which quantifies steepness:

$$n_H=\frac{log81}{log{K}_{90}/K_{10}}\text{,}$$

where K₉₀ is the value of K such that R^ES= − 0.9 and K₁₀ is the value of K such that R^ES= − 0.1. We obtain n _H  = 1.02 when [S] _T  = 10000nM, n _H  = 0.98 when [S] _T  = 100nM and n _H  = 0.73 when [S] _T  = 1nM. Thus, the steepness slightly decreases as the total concentration decreases. Furthermore, we calculate the inhibitor concentration dependence of the complex concentration (Figure 1B). Let the value of the total inhibitor concentration necessary to reach a 50% reduction of concentration (IC₅₀) be a measure of the inhibition. We obtain IC₅₀ = 67.89nM when K = 10nM, IC₅₀ = 59.72nM when K = 100nM and IC₅₀ = 53.20nM when K = 1000nM. We can see that the inhibition becomes more efficient as the equilibrium constant gets larger. Therefore, we conclude that the response coefficient can be a measure of the inhibitory efficiency.

Figure 1

Behaviour of the steady state for a simple model. (A) Response coefficient R as a function of the equilibrium constant K for total substrate concentrations [S] _T  = 100nM, 1000nM and 10000nM. (B) Complex concentration [ES] as a function of the total inhibitor concentration [I] _T for equilibrium constant values K = 10nM, 100nM and 1000nM.

So far, we have discussed behaviours of the response coefficient in a simple model as an example. We would like to generalize this analysis to larger systems which have multiple interactions. From here, we consider models for general classes of complex formation processes: (1) chain reaction systems composed of ordered steps, (2) chain reaction systems composed of unordered multi-branched steps and (3) site-binding reaction systems composed of unordered multi-branched steps.

Analysis of the ordered step models

We first investigate models for chain reaction systems composed of ordered steps. Let A₁, …, A _n be n types of molecules. The reaction scheme of the model is described as

$$\begin{array}{l}{A}_1+A_2\begin{array}{c}\hfill \stackrel{k_{1A}}{\to }\hfill \\ \hfill \underset{k_{1d}}{\leftarrow }\hfill \end{array}{A}_1{A}_2,\\ A_1{A}_2+A_3\begin{array}{c}\hfill \stackrel{k_{2A}}{\to }\hfill \\ \hfill \underset{k_{2d}}{\leftarrow }\hfill \end{array}{A}_1{A}_2{A}_3,\\ A_1{A}_2{A}_3+A_4\begin{array}{c}\hfill \stackrel{k_{3A}}{\to }\hfill \\ \hfill \underset{k_{3d}}{\leftarrow }\hfill \end{array}{A}_1{A}_2{A}_3{A}_4,\\ ⋮\ddots \\ A_1\cdots A_{n-1}+A_n\begin{array}{c}\hfill \stackrel{k_{n-1A}}{\to }\hfill \\ \hfill \underset{k_{n-1d}}{\leftarrow }\hfill \end{array}{A}_1\cdots A_{n-1}{A}_n\text{.}\end{array}$$

We refer to this model as “Model Sn”. In order to provide a basis for the analysis, we begin with the simple case Model S3, which consists of three molecules: A₁, A₂, A₃. By linear approximation under the assumption [A₁] < < [A₂], [A₃], a simple analytical steady state solution can be derived for molecular concentrations

$$\begin{array}{l}\left[A_1\right]=\frac{{\tilde{K}}_1{\tilde{K}}_2}{\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1}{\left[A_1\right]}_T,\\ \left[A_1{A}_2\right]=\frac{{\tilde{K}}_2}{\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1}{\left[A_1\right]}_T,\\ \left[A_1{A}_2{A}_3\right]=\frac{{\tilde{K}}_2}{\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1}{\left[A_1\right]}_T\text{,}\end{array}$$

(2)

where ${\tilde{K}}_1=K_1/{\left[A_2\right]}_T$ and ${\tilde{K}}_2=K_2/{\left[A_3\right]}_T$. The response coefficients for the complex concentrations with respect to K₁ and K₂ are calculated as

$$\begin{array}{l}{R}_1^{123}=\frac{\partial ln\left[A_1{A}_2{A}_3\right]}{\partial ln{\tilde{K}}_1}=-\frac{{\tilde{K}}_2{\tilde{K}}_1}{\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1},\\ R_2^{123}=\frac{\partial ln\left[A_1{A}_2{A}_3\right]}{\partial ln{\tilde{K}}_2}=-\frac{{\tilde{K}}_2\left({\tilde{K}}_1+1\right)}{\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1},\\ R_1^{12}=\frac{\partial ln\left[A_1{A}_2\right]}{\partial ln{\tilde{K}}_1}=-\frac{{\tilde{K}}_2{\tilde{K}}_1}{\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1},\\ R_2^{12}=\frac{\partial ln\left[A_1{A}_2\right]}{\partial ln{\tilde{K}}_2}=\frac{1}{\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1},\\ R_1^1=\frac{\partial ln\left[A_1\right]}{\partial ln{\tilde{K}}_1}=\frac{{\tilde{K}}_2+1}{\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1},\\ R_2^1=\frac{\partial ln\left[A_1\right]}{\partial ln{\tilde{K}}_2}=\frac{1}{\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1}\text{,}\end{array}$$

(3)

where R _i ¹²³, R _i ¹² and R _i ¹(i = 1, 2) denote the response coefficients for [A₁A₂A₃], [A₁A₂] and [A₁] with respect to K _i , respectively. As can be seen from (3), the following inequalities between the response coefficients hold: R₁¹²³ > R₂¹²³, R₂¹² > R₁¹² and R₁¹ > R₂¹. The inequality R₁¹²³ > R₂¹²³ shows that inhibitory efficiency of the A₃ binding reaction is larger than that of the A₂binding reaction. The sign of R₁¹² is negative, while that of R₂¹² is positive. This means that the complex concentration [A₁A₂] increases by the inhibition of the A₂ binding reaction.

In the nonlinear regime of Model S3, we cannot solve for the steady state solution explicitly. Thus, we instead investigate whether the same relations hold by numerical simulations. In the simulations, we set the total concentrations as [A _i ] _T  = 100nM(i = 1, 2, 3). From the simulation results, we can see the following: R₁¹²³ depends on the values of K₁and K₂ (Figure 2A). The inequality R₁¹²³ > R₂¹²³ holds for a wide range of K₁ values (Figure 2B). The sign of R₁¹² is negative, while that of R₂¹² is positive (Figure 2C). R₁¹ and R₂¹ are both positive and satisfy R₁¹ > R₂¹ (Figure 2D). Therefore, we have shown that the parameter-independent inequalities hold both in the linear and the nonlinear regimes of Model S3.

Figure 2

Behaviours of the response coefficients for Model S3. (A) Response coefficient R₁¹²³ as a function of K₁ for K₂ = 100nM, 1 nM and 0.01 nM. (B) Response coefficient R _i ¹²³ (i = 1, 2) as a function of K₁. Blue and red lines indicate R₁¹²³ and R₂¹²³, respectively. (C) Response coefficient R _i ¹²(i = 1, 2) as a function of K₁. Blue and red lines indicate R₁¹² and R₂¹², respectively. (D) Response coefficient R _i ¹(i = 1, 2) as a function of K₁. Blue and red lines indicate R₁¹ and R₂¹, respectively. In (B)–(D), we set K₂ = 100nM.

To investigate whether similar relations hold in larger systems, we next consider Model S4. To linearize the model, we assume [A₁] _T  < < [A₂] _T , [A₃] _T , [A₄] _T . As in Model S3, the steady state solution can be derived for molecular concentrations

$$\begin{array}{l}\left[A_1\right]=\frac{{\tilde{K}}_1{\tilde{K}}_2{\tilde{K}}_3}{\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right){\tilde{K}}_3+1}{\left[A_1\right]}_T,\\ \left[A_1{A}_2\right]=\frac{{\tilde{K}}_2{\tilde{K}}_3}{\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right){\tilde{K}}_3+1}{\left[A_1\right]}_T,\\ \left[A_1{A}_2{A}_3\right]=\frac{{\tilde{K}}_3}{\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right){\tilde{K}}_3+1}{\left[A_1\right]}_T,\\ \left[A_1{A}_2{A}_3{A}_4\right]=\frac{1}{\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right){\tilde{K}}_3+1}{\left[A_1\right]}_T\text{.}\end{array}$$

(4)

The response coefficients are calculated as

$$\begin{array}{l}{R}_1^{1234}=-\frac{{\tilde{K}}_3{\tilde{K}}_2{\tilde{K}}_1}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1}\text{,}{R}_2^{1234}=-\frac{{\tilde{K}}_3{\tilde{K}}_2\left({\tilde{K}}_1+1\right)}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1}\text{,}\\ R_3^{1234}=-\frac{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1},R_1^{123}=-\frac{{\tilde{K}}_3{\tilde{K}}_2{\tilde{K}}_1}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1}\text{,}\\ R_2^{123}=-\frac{{\tilde{K}}_3{\tilde{K}}_2\left({\tilde{K}}_1+1\right)}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1}\text{,}{R}_3^{123}=\frac{1}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1},\\ R_1^{12}=-\frac{{\tilde{K}}_3{\tilde{K}}_2{\tilde{K}}_1}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1}\text{,}{R}_2^{12}=\frac{{\tilde{K}}_3+1}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1}\text{,}\\ R_3^{12}=\frac{1}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1},R_1^1=\frac{\left({\tilde{K}}_2+1\right){\tilde{K}}_3+1}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1}\text{,}\\ R_2^{12}=\frac{{\tilde{K}}_3+1}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1}\text{,}{R}_3^{12}=\frac{1}{{\tilde{K}}_3\left(\left({\tilde{K}}_1+1\right){\tilde{K}}_2+1\right)+1}\text{,}\end{array}$$

(5)

where $R_i^{1234}=\frac{\partial ln\left[A_1{A}_2{A}_3{A}_4\right]}{\partial ln{\tilde{K}}_i}$, $R_i^{123}=\frac{\partial ln\left[A_1{A}_2{A}_3\right]}{\partial ln{\tilde{K}}_i}$, $R_i^{12}=\frac{\partial ln\left[A_1{A}_2\right]}{\partial ln{\tilde{K}}_i}$ and $R_i^1=\frac{\partial ln\left[A_1\right]}{\partial ln{\tilde{K}}_i}(i=1,2,3)$. From the above expression (5), we can see that the inequalities R₁¹²³⁴ > R₂¹²³⁴ > R₃¹²³⁴, R₃¹²³ > R₁¹²³ > R₂¹²³, R₂¹² > R₃¹² > R₁¹² and R₁¹ > R₂¹ > R₃¹ hold. In the simulations, we set total concentrations [A _i ] _T  = 100nM(i = 1, 2, 3, 4). We can see that the inequalities R₁¹²³⁴ > R₂¹²³⁴ > R₃¹²³⁴, R₃¹²³ > R₁¹²³ > R₂¹²³, R₂¹² > R₃¹² > R₁¹² and R₁¹ > R₂¹ > R₃¹ hold for a wide range of K₁ values (Figure 3). The sign of R₁¹²³ is positive, while those of R₂¹²³ and R₃¹²³ are negative (Figure 3B). R₁¹² and R₂¹² are positive, while R₃¹² is negative (Figure 3C). R₁¹, R₂¹ and R₃¹ are positive (Figure 3D). These results suggest that the inequalities relating the response coefficients represent basic properties of the ordered step models.

Figure 3

Behaviours of the response coefficients for Model S4. (A) Response coefficient R _i ¹²³⁴ (i = 1,2,3) as a function of K₁. Blue, red and green lines indicate R₁¹²³⁴, R₂¹²³⁴ and R₃¹²³⁴, respectively. (B) Response coefficient R _i ¹²³(i = 1, 2, 3) as a function of K₁. Blue, red and green lines indicate R₁¹²³, R₂¹²³ and R₃¹²³, respectively. (C) Response coefficient R _i ¹²(i = 1, 2, 3) as a function of K₁. Blue, red and green lines indicate R₁¹², R₂¹² and R₃¹², respectively. (D) Response coefficient R _i ¹(i = 1, 2, 3) as a function of K₁. Blue, red and green lines indicate R₁¹, R₂¹ and R₃¹, respectively. For all, we set K₂ = K₃ = 100nM.

We can generalize these results to the n chain reaction model: Model Sn. Assuming [A₁] _T  < < [A₂] _T , …, [A _n ] _T , the steady state solution can be obtained as

$$\begin{array}{l}\left[A_1\right]=\frac{{\tilde{K}}_{n-1}\cdots {\tilde{K}}_1}{D}{\left[A_1\right]}_T,\\ ⋮\\ [A_1\cdots A_i]=\frac{{\tilde{K}}_{n-1}\cdots {\tilde{K}}_i}{D}{\left[A_1\right]}_T\\ ⋮\\ [A_1{A}_2\cdots A_n]=\frac{1}{D}{\left[A_1\right]}_T\text{,}\end{array}$$

(6)

where $D=1+{\tilde{K}}_{n-1}\left(1+{\tilde{K}}_{n-2}\left(1+\cdots \left(1+{\tilde{K}}_1\right)\right)\right)$. The response coefficients are calculated as

$$\begin{array}{l}{R}_{n-1}^{1\cdots i}=\frac{1}{D},\\ ⋮\\ R_{i+1}^{1\cdots i}=\frac{1+{\tilde{K}}_{n-1}\left(1+\cdots \left(1+{\tilde{K}}_{i+2}\right)\right)}{D},\\ R_i^{1\cdots i}=-\frac{{\tilde{K}}_{n-1}\cdots {\tilde{K}}_i\left(1+{\tilde{K}}_{i-1}\left(1+\cdots \left(1+{\tilde{K}}_1\right)\right)\right)}{D},\\ ⋮\\ R_1^{1\cdots i}=-\frac{{\tilde{K}}_{n-1}\cdots {\tilde{K}}_1}{D}\text{.}\end{array}$$

(7)

where $R_i^j=\frac{\partial ln[A_1\cdots A_j]}{\partial ln{\tilde{K}}_i}$. Therefore, the inequalities are as follows:

$$\begin{array}{l}{R}_1^{1\cdots n}>R_2^{1\cdots n}>\cdots >R_{n-1}^{1\cdots n},\\ ⋮\\ R_{i+1}^{1\cdots i}>R_{i+2}^{1\cdots i}>\cdots >R_{n-1}^{1\cdots i}>R_1^{1\cdots i}>\cdots >R_i^{1\cdots i},\\ ⋮\\ R_1^1>R_2^1>\cdots >R_{n-1}^1\text{.}\end{array}$$

Thus, we have shown that the parameter-independent inequalities hold in Model Sn. The inequalities indicate that the reactions can be sorted in order of the inhibitory efficiency, which is independent of the values of the equilibrium constants.

Analysis of the unordered multi-branched step models – chain reaction systems

Here, we consider chain reaction systems with unordered multi-branched steps. In the unordered multi-branched step models, we assume that there are no cooperative reactions, i.e. the affinity of any reaction is independent of whether other sites are occupied. Let A₁, …, A _n be n types of molecules. All possible complexes in the general chain reaction system are listed as follows:

$$\begin{array}{l}{A}_1{A}_2{A}_3\cdots A_{n-1}{A}_n\\ A_1{A}_2{A}_2{A}_3\cdots \cdots \cdots A_{n-1}{A}_n\\ A_1{A}_2{A}_3\cdots \cdots A_{n-2}{A}_{n-1}{A}_n\\ \ddots ⋰\\ A_1{A}_2\cdots A_n\end{array}$$

We refer to this model as “Model Cn”. The number of complexes is N _S  = n(n + 1)/2 and the number of reactions is N _I  = n(n − 1)(n + 1)/2. As an example, let us first consider the n = 3 case. All possible complexes can be listed as

$$\begin{array}{l}{A}_1{A}_2{A}_3\\ A_1{A}_2{A}_2{A}_3\\ A_1{A}_2{A}_3\end{array}$$

In order to deal with the model analytically, we linearize the system. By linear approximation under the assumption [A₂] _T  < < [A₁]_T,[A₃] _T , the steady state solution for molecular concentrations is written simply as

$$\begin{array}{l}\left[A_2\right]=\frac{{\tilde{K}}_1{\tilde{K}}_2}{\left({\tilde{K}}_1+1\right)\left({\tilde{K}}_2+1\right)}{\left[A_2\right]}_T,\\ \left[A_1{A}_2\right]=\frac{{\tilde{K}}_2}{\left({\tilde{K}}_1+1\right)\left({\tilde{K}}_2+1\right)}{\left[A_2\right]}_T,\\ \left[A_2{A}_3\right]=\frac{{\tilde{K}}_1}{\left({\tilde{K}}_1+1\right)\left({\tilde{K}}_2+1\right)}{\left[A_2\right]}_T,\\ \left[A_1{A}_2{A}_3\right]=\frac{1}{\left({\tilde{K}}_1+1\right)\left({\tilde{K}}_2+1\right)}{\left[A_2\right]}_T\text{,}\end{array}$$

(8)

where ${\tilde{K}}_1=K_1/{\left[A_1\right]}_T$, ${\tilde{K}}_2=K_2/{\left[A_3\right]}_T$. The response coefficients for molecular concentrations with respect to the equilibrium constants are given by

$$\begin{array}{l}{R}_1^{123}=-\frac{{\tilde{K}}_1}{{\tilde{K}}_1+1}\text{,}{R}_2^{123}=-\frac{{\tilde{K}}_2}{{\tilde{K}}_2+1},\\ R_1^{12}=-\frac{{\tilde{K}}_1}{{\tilde{K}}_1+1}\text{,}{R}_2^{12}=\frac{1}{{\tilde{K}}_2+1},\\ R_1^{23}=\frac{1}{{\tilde{K}}_1+1}\text{,}{R}_2^{23}=-\frac{{\tilde{K}}_2}{{\tilde{K}}_2+1},\\ R_1^2=\frac{1}{{\tilde{K}}_1+1}\text{,}{R}_2^2=\frac{1}{{\tilde{K}}_2+1}\text{,}\end{array}$$

(9)

where $R_i^{123}=\frac{\partial ln\left[A_1{A}_2{A}_3\right]}{\partial ln{\tilde{K}}_i}$, $R_i^{12}=\frac{\partial ln\left[A_1{A}_2\right]}{\partial ln{\tilde{K}}_i}$, $R_i^{23}=\frac{\partial ln\left[A_2{A}_3\right]}{\partial ln{\tilde{K}}_i}$ and $R_i^2=\frac{\partial ln\left[A_2\right]}{\partial ln{\tilde{K}}_i}(i=1,2)$. An important observation is that the response coefficients with respect to ${\tilde{K}}_i$ only depend on${\tilde{K}}_i$, i.e. $R_i=R_i\left({\tilde{K}}_i\right)$. We call this property independence of the response coefficients. This property comes from the assumption that there are no cooperative reactions, because, as shown previously, the ordered step models do not have this property. Numerical simulation results also show independence of the response coefficients in the nonlinear regime of the model (Figure 4). In the simulations, we set [A _i ] _T  = 100nM(i = 1, 2, 3). The response coefficient R₁¹²³ does not depend on K₂ for a wide range of K₁ values (Figure 4A). The sign of the response coefficient R₂¹²is positive, which means [A₁A₂] increases by the inhibition of K₂ (Figure 4B). The independence of R₂¹² in terms of K₂also holds for a wide range of K₁ values (Figure 4B). The independence of response coefficients means that inhibitory efficiency is determined by the values of the equilibrium constants.

Figure 4

Behaviours of the response coefficients for Model C3. (A) Response coefficient R₁¹²³ as a function of K₁ for K₂ = 100nM, 1nM and 0.01nM. (B) Response coefficient R₂¹² as a function of K₂ for K₁ = 100nM, 1nM and 0.01nM.

We next explore Model C4. All possible complexes are listed as

$$\begin{array}{l}{A}_1{A}_2{A}_3{A}_4\\ A_1{A}_2{A}_2{A}_3{A}_3{A}_4\\ A_1{A}_2{A}_3{A}_2{A}_3{A}_4\\ A_1{A}_2{A}_3{A}_4\end{array}$$

We investigate this model only by means of numerical simulations. In the simulations, we set [A _i ] _T  = 100nM(i = 1, 2, 3, 4). The response coefficient R₁¹²³⁴ does not depend on K₂ and K₃ for a wide range of K₁ values (Figure 5A). The independence of R₃¹²³ in terms of K₁ and K₂ also holds for a wide range of K₃ values (Figure 5B).

Figure 5

Behaviours of the response coefficients for Model C4. (A) Response coefficient R₁¹²³⁴ as a function of K₁ for different K₂ and K₃ values. Blue (control) line indicates K₂ = K₃ = 100nM, red line indicates K₂ = 1nM, K₃ = 100nM, and green line indicates K₂ = 100nM, K₃ = 1nM. (B) Response coefficient R₃¹²³ as a function of K₃ for different K₁ and K₂ values. Blue (control) line indicates K₁ = K₂ = 100nM, red line indicates K₁ = 1nM, K₂ = 100nM, and green line indicates K₁ = 100nM, K₂ = 1nM.

Analysis of the unordered multi-branched step models – site-binding systems

Here, we apply control analysis to the models in which molecule A has n binding sites for the binding of B _i (i = 1, ⋯, n). Schematically, all possible complexes are listed as

$$\begin{array}{l}A{B}_1{B}_2\cdots B_n\\ {\text{AB}}_1{\text{AB}}_2\cdots {\text{AB}}_n\\ {\text{AB}}_1{B}_2{\text{AB}}_1{B}_3\cdots {\text{AB}}_{n-1}{B}_n\\ \ddots ⋮\\ {\text{AB}}_1\cdots B_{n-1}\cdots {\text{AB}}_2\cdots B_n\\ {\text{AB}}_1\cdots B_n\end{array}$$

We refer to this model as “Model SBn”. We assume that each equilibrium constant K _i for binding the site of B_ito each site of A is independent of whether any of the other sites of A are occupied. By linear approximation under the assumption [A] _T  < < [B₁] _T , ⋯, [B _n ] _T , the steady state solution for molecular concentrations is given by

$$\begin{array}{l}\left[A\right]=\left(\prod _{j=1}^n{K}_j\right)A_{1\cdots n},\\ \left[A{B}_i\right]=\left(\prod _{j\ne i}^n{K}_j\right)A_{1\cdots n}\text{,}i=1,\cdots ,n\\ ⋮\\ [A{B}_{i1}\cdots B_{\mathit{ik}}]=\left(\prod _{j\ne i1,\cdots ,ik}^n{K}_j\right)A_{1\cdots n}\text{,}i1,\cdots ,ik(i1\ne \cdots \ne ik)=1,\cdots ,n\\ ⋮\\ [A{B}_1\cdots B_n]=A_{1\cdots n}\text{,}\end{array}$$

(10)

where $A_{1\cdots n}=\left(\prod _{j=1}^n(1+K_j)\right)$. The response coefficients for [AB_{i 1}⋯B _ik ] with respect to K _j are given by

$$R_j^{i1\cdots ik}=\begin{array}{c}\hfill -\frac{{\tilde{K}}_j}{1+{\tilde{K}}_j}ifj=i1,\cdots ,ik(i1\ne \cdots \ne ik),\hfill \\ \hfill \frac{1}{1+{\tilde{K}}_j}ifj\ne i1,\cdots ,ik(i1\ne \cdots \ne ik)\text{,}\hfill \end{array}$$

(11)

where $R_j^{i1\cdots ik}=\frac{\partial ln[A_{i1}\cdots A_{\mathit{ik}}]}{\partial ln{K}_j}$. Thus, the independence of the response coefficients also holds in this model. We further explore the nonlinear regime of the model by using numerical simulations. We consider Model SB3 as an example. In the simulations, we set [A] _T  = [B₁] _T  = ⋯ = [B _n ] _T  = 100nM. We can see that the independence of R₁¹²³ holds in a wide range of K₁ (Figure 6). Therefore, together with the results from the chain reaction systems, the independence of the response coefficients is an intrinsic property of unordered multi-branched step models under the assumption that there are no cooperative reactions.

Figure 6

Behaviour of the response coefficient for Model SB3. Response coefficient R₁¹²³ as a function of K₁ for K₂ = 100nM, 1nM and 0.01nM. Here, we set K₃ = 100nM

An application - MT1-MMP/TIMP2/MMP2 complex formation model

As an application of our analysis, we here investigate the MT1-MMP/TIMP2/MMP2 complex formation model [16]. In this model, MT1-MMP is dimerized, bound to TIMP2 and forms a quadruple complex by binding proMMP2, MT1-MT1T2M2, which contributes to activation of proMMP2. The model has all possible pathways that could contribute to the formation and dissociation of complexes. All possible complexes formed by the mentioned biochemical reactions are listed below:

$$\begin{array}{l}\text{MT}1\text{T}2\text{M}2\\ \text{MT1-MT}1\text{MT1T}2\text{T2M}2\\ \text{MT1-MT1T}2\text{MT1T2M}2\\ \text{MT1T2-MT1T}2\text{MT1-MT1T2M}2\\ \text{MT1T2-MT1T2M}2\\ \text{MT1T2M2-MT1T2M}2\end{array}$$

Rate constants of the reactions are summarized in Table 1. Although a simpler model was proposed previously in [18], we use the model described above. As reported in [16], taking into consideration the transient dynamics of the MMP2 activation, the model we use here is more appropriate than the model proposed in [18].

Table 1 Experimentally derived parameter values for the MMP model

Here, we try to identify the most efficient reaction to consider in selecting inhibitors by putting the experimentally derived parameter values into the model. We focus on the behaviour of the concentration of the quadruple complex MT1-MT1T2M2, which is a binding state of a free MT1-MMP and the ternary complex MT1T2M2, because only this complex contributes to activation of proMMP2. The total TIMP2 concentration dependence of the complex MT1-MT1T2M2 is shown (Figure 7A). The concentration of the quadruple complex has the maximum value at a TIMP-2 concentration of 50 nM (Figure 7B control). This property does not depend on MT 1] _T and M 2] _T . This means that the quadruple complex which contributes to the proMMP2 activation is suppressed at high levels of TIMP2 and attains its maximum only at intermediate TIMP2 concentrations [16–18]. A concentration map in the MT 1] _T and T 2] _T plane is depicted in Figure 7B. These simulation results indicate that the balance between TIMP2 and MT1-MMP expression is a critical determinant of MMP2 activation. The response coefficients for MT1-MT1T2M2 with respect to the equilibrium constants were also calculated (Figure 7C). This model assumes that there are no cooperative reactions, as in the previous examples, and thus the independence of the response coefficients holds. To compare the three curves R_{MT 1MT 1}, R_{MT 1T 2} and R_{T 2M 2} quantitatively, we calculate the steepness. The Hill coefficients are n _H  = 0.74 for R_{MT 1MT 1}, n _H  = 1.15 for R_{MT 1T 2} and n _H  = 0.95 for R_{T 2M 2}. Thus, their steepnesses are slightly different. By substituting experimentally derived values of the equilibrium constant (Table 2) into the response coefficients, we obtain R_{MT 1MT 1} = − 0.095(K_{MT 1MT 1} = 5nM), R_{MT 1T 2} = − 0.0024(K_{MT 1T 2} = 0.548nM) and R_{T 2M 2} = − 0.15(K_{T 2M 2} = 33.5714nM). Therefore, we conclude that the most effective inhibition is the TIMP2 and MMP2 binding interaction.

Figure 7

Behaviour of the steady state of the MMP model. (A) Complex concentration [MT 1MT 1T 2M 2] as a function of total TIMP2 concentration for different total MT1-MMP and MMP2 concentrations. Blue line indicates [MT 1] _T  = 100nM, [M 2] _T  = 100nM. Red line indicates [MT 1] _T  = 50nM, [M 2] _T  = 100nM. Green line indicates [MT 1] _T  = 10nM, [M 2] _T  = 100nM. Purple line indicates [MT 1] _T  = 100nM, [M 2] _T  = 50nM. Cyan line indicates [MT 1] _T  = 100nM, [M 2] _T  = 10nM. (B) Concentration map of the complex concentration [MT 1MT 1T 2M 2] as a function of total MT1-MMP and TIMP2 concentrations. Here, we set [M 2] _T  = 100nM. (C) The response coefficients of the MMP model for the three interactions: the MT1-MMP dimerization, the MT1-MMP and TIMP2 binding, and the TIMP2 and MMP2 binding. The blue line shows the response coefficient R_{MT 1MT 1} = ∂ ln [MT 1MT 1T 2M 2]/ ∂ ln K_{MT 1MT 1} as a function of K_{MT 1MT 1}, the red line shows the response coefficient R_{MT 1T 2} = ∂ ln [MT 1MT 1T 2M 2]/ ∂ ln K_{MT 1T 2} as a function of K_{MT 1T 2} and the green line shows the response coefficient R_{T 2M 2} = ∂ ln [MT 1MT 1T 2M 2]/ ∂ ln K_{T 2M 2} as a function of K_{T 2M 2}. Here, we set [MT 1] _T  = 100nM, [T 2] _T  = 50nM and [M 2] _T  = 100nM.

Table 2 Summary of equilibrium constants

Several broad-spectrum MMPIs function by strongly chelating the Zn ion that lies in the MMP active site. Similarities in active sites of MMPs pose obstacles to the design of specific inhibitors [3]. Thus, a new approach to the identification of new drug targets is important. Here, focusing on the MMP2 activation process, we were able to determine that the TIMP2-MMP2 is the most effective interaction. It is reported that TIMP2 interacts with MMP2 through the C-terminal domain of the enzyme that is distinct from the active site [21, 22]. Therefore, our result identifies a new drug target in the process of MMP2 activation. Development of low molecular weight compounds capable of effectively and specifically inhibiting the TIMP2 and MMP2 binding interaction will be the subject of future research. Our result can be validated using cell culture systems.

Conclusions

In this paper, our aim is to quantify the response of a system to the addition of inhibitors and to classify their interactions in order of their inhibitory efficiency. In order to analyse the response systematically, we used control analysis. Using the response coefficients, we revealed that the parameter-independent inequalities between the response coefficients hold in the ordered step models. For the unordered multi-branched step models, we showed that independence of the response coefficients with respect to the equilibrium constants holds. These results indicate that the inhibitory efficiency depends on the topology of the pathway networks. We applied our analysis to a complex formation model describing the formation of complexes of MMP2 and MT1-MMP in the presence of TIMP2 [16]. In the complex formation process between these molecules, there are three interactions, i.e. the MT1-MMP dimerization, the MT1-MMP and TIMP2 binding reaction and the TIMP2 and MMP2 binding reaction. We tried to identify the most efficient interaction to consider in selecting the inhibitors by putting the experimentally derived parameter values into the model. The novel finding of the analysis is that the inhibition of the TIMP2 and MMP2 binding interaction is the most efficient method for suppressing the quadruple complex MT1-MT1T2M2, which contributes to the MMP2 activation. This result identifies a new drug target in the process of MMP2 activation.

Our method can also be applied to other models of complex formation processes. However, there are some weaknesses in the analysis presented here. Throughout our treatment, we have considered only the steady state in a well-stirred environment. For the case that the substrate is a non-diffusible molecule such as ECM, the well-stirred assumption is not valid. Thus, in this case, we should consider a three dimensional compartment model in which the extracellular space is divided into small compartments [16, 23]. This will allow us to simulate a reaction–diffusion system. Furthermore, we analysed models in closed systems, but intra- and intercellular transport of molecules should play important roles in biochemical reactions occurring in cells. Thus, spatio-temporal dynamics associated with a mechanism such as a positive feedback loop [24] should be considered.

Methods

Numerical computation scheme

We employed the fourth-order Runge–Kutta method to solve systems of ordinary differential equations numerically. In all simulations, the time step was taken as dt = 0.001 sec and time evolution was performed up to the time T = 100000 sec. In the calculation of the response coefficients, the small fractional change of the equilibrium constant K was δK/K = 0.01.