A Biochemical Model of Matrix Metalloproteinase 9 Activation and Inhibition*

Matrix metalloproteinases (MMPs) are a class of extracellular and membrane-bound proteases involved in an array of physiological processes, including angiogenesis. We present a detailed computational model of MMP9 activation and inhibition. Our model is validated to existing biochemical experimental data. We determine kinetic rate constants for the processes of MMP9 activation by MMP3, MMP10, MMP13, and trypsin; inhibition by the tissue inhibitors of metalloproteinases (TIMPs) 1 and 2; and MMP9 deactivation. This computational approach allows us to investigate discrepancies in our understanding of the interaction of MMP9 with TIMP1. Specifically, we find that inhibition due to a single binding event cannot describe MMP9 inhibition by TIMP1. Temporally accurate biphasic inhibition requires either an additional isomerization step or a second lower affinity isoform of MMP9. We also theoretically characterize the MMP3/TIMP2/pro-MMP9 and MMP3/TIMP1/pro-MMP9 systems. We speculate that these systems differ significantly in their time scales of activation and inhibition such that MMP9 is able to temporarily overshoot its final equilibrium value in the latter. Our numerical simulations suggest that the ability of pro-MMP9 to complex TIMP1 increases this overshoot. In all, our analysis serves as a summary of existing kinetic data for MMP9 and a foundation for future models utilizing MMP9 or other MMPs under physiologically well defined microenvironments.

Matrix metalloproteinase 9 (MMP9) 2 is a zinc-dependent endopeptidase that participates in a variety of physiological and biochemical processes (1,2). In angiogenesis, MMP9 has been shown to play a part in the degradation of basement membranes through its type IV collagenase activity (3), creation of bioactive molecules such as angiostatin by cleavage of plasminogen (4), and the release of growth factors such as vascular endothelial growth factor-A sequestered on extracellular matrix heparin sulfate proteoglycans (5)(6)(7). To quantify the existing biochemical data and help understand properties of MMP activation cascades, we have constructed a biochemically accurate numerical model of the post-secretional regulation of MMP9. Specifically, we have used existing experimental data to kinetically characterize mechanistic models of the MMP9 reaction network (e.g. Fig. 1). MMP9 enters this model (post-secretion) as a glycosylated latent enzyme (pro-MMP9) at a molecular mass of 92 kDa. From experimental studies, MMP9 is known to undergo activation, inhibition by tissue inhibitors of metalloproteinases (TIMPs), and possibly post-activation processing, including C-terminal processing and deactivation. In this approach, we have not studied the proteolytic action elicited by MMP9. A vast quantity of literature exists to create detailed models of activation and inhibition; however, C-terminal processing and MMP9 deactivation are less studied or accepted.
In general, MMP latency is maintained by an N-terminal prodomain that masks the Zn 2ϩ -dependent active site. The structure of the active molecule and the specific rate of activation depend on the exact mode of activation, but in all cases they require conformational unwinding in the catalytic domain. By far, the most studied mode of MMP9 activation is enzyme proteolysis of the pro-MMP9 prodomain. Confirmed in vitro activators of MMP9 include MMP2 (8), MMP3 (9), MMP7 (10), MMP10 (11), MMP13 (12), and the serine protease trypsin (13). MMP3 is generally regarded as the most potent activator of MMP9. Here, MMP9 is processed via an intermediate of 86 kDa (processing at Glu 59 -Met 60 ), which relaxes to expose the second processing site (Arg 106 -Phe 107 ) leading to the well established 82-kDa isoform. Evidence also supports nonproteolytic modes of activation, e.g. via conformational changes contingent on substrate binding (14). Post-activation processing for MMP9 can further involve removal of the C-terminal hemopexin-like domain, leading to an ϳ65-kDa species (MMP9⌬C), or can involve direct cleavage in the active domain producing an ϳ50 -60-kDa inactive species (15). Thus, MMP9 is typically processed from 92 to 86 kDa and 82, 65, and ϳ50 -60 kDa. A distinct feature of MMP9 is its ability to form proenzyme complexes with TIMP1 (pro-MMP9⅐TIMP1) and to form dimers with other MMPs, including the MMP9⅐MMP1 heterodimer (16) and MMP9⅐MMP9 homodimer (17).
Once activated, MMP9 is subject to inhibition by TIMPs resulting from binding of the N-terminal domain of a TIMP to the active site cleft in MMP9. Our understanding of MMP9 inhibition is incomplete. First, MMP9 is found as active isoforms at both 82 and 65 kDa, but the 65-kDa isoform lacks the hemopexin-like domain, which is generally recognized to increase the rate of MMP9 inhibition by TIMP1 (18). We do not know what role enzymatic cleavage of the C-terminal domain of MMP9 plays in MMP9 inhibition. In addition, although MMP TIMP pairs such as MMP9 and TIMP2 or MMP2 and TIMP1 exhibit monophasic inhibition kinetics, MMP TIMP combinations capable of forming pro-MMP⅐TIMP complexes (e.g. pro-MMP9⅐TIMP1 and pro-MMP2⅐TIMP2) seem to exhibit biphasic inhibition kinetics (18,19). An important problem in our current understanding of MMP9 inhibition by TIMP1 is the large discrepancy in the reported K i values between various groups. Assuming simple bimolecular kinetics, O'Connell et al. (18) reported values of K i Ͻ50 pM, with an association rate constant (k on ) of 1.11 ϫ 10 7 M Ϫ1 s Ϫ1 . In contrast, Olson et al. (19) reported lower affinity kinetics with nanomolar K i ϭ 8.5 nM and k on equal to 2.48 ϫ 10 5 M Ϫ1 s Ϫ1 . Furthermore, resonance biosensor analyses between MMPs and TIMPs consistently indicate binding kinetics with K i in the low and intermediate nanomolar range across a wide range of MMP and TIMP combinations (19,20). An important goal of this study is to examine the existing relationships and use numerical simulations to determine models capable of accurately reproducing observed behavior.
In the physiological setting, TIMPs play a major role in inhibiting MMPs, and thus pro-MMP9 activation will never truly occur in isolation of an inhibitor. Experiments indicate that the rate of pro-MMP9 activation is reduced in the presence of TIMPs (15,21) because of either TIMP1 inhibition of the activators or TIMP1 stabilization of pro-MMP9 by forming the pro-MMP9⅐TIMP1 complex. TIMP1 can bind pro-MMP9 through its C-terminal domain, leaving its active site available to inhibit other MMPs, e.g. MMP3. This is responsible for the importance of the MMP3 to TIMP1 ratio on pro-MMP9⅐TIMP1 complex activation.
In all, our computational approach allowed us to test various hypotheses regarding MMPs and TIMPs against previously observed behavior. Using similar methodology, we previously confirmed the nonmonotonic dependence of MMP2 activation on MT1-MMP and TIMP2 ratios in the MT1-MMP/TIMP2/ MMP2 system and showed the importance of membranebound MMPs over diffusible MMPs in type 1 collagen degradation (22).

MATERIALS AND METHODS
Based on the mechanisms proposed in the literature, we formulated a network of reactions in which MMP9 and its activators and inhibitors are likely to participate (Figs. 1 and 3). To form a numerical model of MMP9 behavior, these reactions were converted to a system of ordinary differential equations by applying mass action kinetics. These equations express the rate of change in concentration of each species based on the production of each species or consumption rates because of the listed reactions. Specifying the initial concentrations of all species allows calculation of the concentration of each species at a later time. For the simulated reactions, we assume well stirred and cell-free conditions allowing us to neglect spatial gradients of the species.
Next, kinetic parameters for these reactions were either obtained from values calculated by other groups or were independently determined using data from the experimental literature (we did not carry out our own biochemical experiments).
To determine kinetic parameters, we employed a grid search algorithm to minimize normalized error (1 Ϫ R 2 ϭ sum of squared residuals normalized to the total sum of squares) between the simulated results of our model and the experimental data (usually of time courses or dose-response curves).
The specific reactions and equations, and the methodology used in estimating unknown kinetic parameters for these reactions, are detailed in the Supplemental Material. The resulting system of ordinary differential equations was solved with the Matlab software package (Mathworks, Inc.) using their ode23s stiff ordinary differential equation solver. Numerical precision was double floating point, and the time step size was automatically selected to have a relative error tolerance less than 10 Ϫ3 and an absolute error tolerance less than 10 Ϫ15 .

RESULTS AND DISCUSSION
Estimates for Pro-MMP9 Activation in the Absence of TIMPs-Pro-MMP9 activation is a sequential processing event involving the conversion of the latent 92-kDa form to an 82-kDa species by an 86-kDa inactive intermediate (9). We treated activation as a single step event leading from a 92-to an 82-kDa form. This assumption is valid because the 92-to 86-kDa conversion usually occurs fairly rapidly (several seconds) compared with activation to the 82-kDa species, which takes hours (18).
First, we numerically estimated pro-MMP9 activation by MMP3 ( Fig. 2A) using data from an experiment by O'Connell et al. (18). For the case of pro-MMP9 activation by MMP3, a Michaelis-Menten activation constant k act was previously determined to be 0.0019 s Ϫ1 (23). By performing a numerical sensitivity analysis, we estimated different k on and k off values that fit the experimental time course to within experimental uncertainty as shown by the contour plot in Fig. 2B. These range from the relatively slow binding kinetics at k on ϭ 6 ϫ 10 3 M Ϫ1 s Ϫ1 and k off ϭ 0 s Ϫ1 to relatively rapid kinetics at k on ϭ 6 ϫ 10 5 M Ϫ1 s Ϫ1 and k off ϭ 0.1 s Ϫ1 . The absolute minimum for k on assuming that k act is infinity is 2.5 ϫ 10 3 M Ϫ1 s Ϫ1 . With k act ϭ 0.0019 s Ϫ1 , we could restrict k on to being at least 6 ϫ 10 3 M Ϫ1 s Ϫ1 . We found that our later simulations were not very sensitive to different parameter combinations; thus we adopted one set of parameters for the rest of this study chosen as the median values of the estimated binding kinetics: k on ϭ 1 ϫ 10 4 M Ϫ1 s Ϫ1 , k off ϭ 1.0 ϫ 10 Ϫ3 s Ϫ1 , and k act ϭ 0.0019 s Ϫ1 (k act /K m ϭ 8840 M Ϫ1 s Ϫ1 ). The entire family of optimal parameters leads to effective k act /K m in the range of 6 ϫ 10 3 M Ϫ1 s Ϫ1 to 1.1 ϫ 10 4 M Ϫ1 s Ϫ1 , which is almost 10 times lower than that determined previously by Olson et al. (23). We also carried out similar analyses for MMP10 (Fig. 2C), MMP13, and trypsin. The results are presented in Table 1. Our calculations indicate that MMP13 and MMP3 are equally competent MMP9 activators, followed by trypsin, then MMP10. To show this, we simulated a test experiment where 200 nM pro-MMP9 was separately combined with 30 nM of each activator (Fig. 2D).

MMP9 Interaction with TIMP1 Cannot Be Described by a Single
Step Binding Inhibition Reaction-To study the interaction of MMP9 and TIMP1, we numerically tested several previously suggested reaction mechanisms against the progress curve data of Olson et al. (19). These progress curves (Fig. 3) followed the onset of inhibition between free MMP9 and both TIMP1 and TIMP2 and the resumption of MMP9 activity following the dissociation of the MMP9⅐TIMP complex.
The experimental data show a clear difference between the inhibitory behavior of TIMP1 and TIMP2 for MMP9. For example, similar increases of TIMP1 or TIMP2 elicit a larger loss of MMP9 activity within the first 100 s and a smaller loss of activity past 100 s for TIMP1. This indicates that TIMP1 causes a biphasic pattern of inhibition with an initial phase of relatively rapid inhibition followed by a phase of slower inhibition.
The four mechanisms we tested include a single step binding inhibition (scheme I), a single step binding inhibition followed by isomerization into a second inhibited complex (scheme II), a noninhibitory initial binding followed by isomerization to an inhibited complex (scheme III), and a two-isoform MMP9 model with each MMP9 isoform having different binding properties (scheme IV). The reaction mechanisms for these schemes are displayed in Fig. 3. For each case, we present the predicted progress curves for the inhibition of MMP9 by TIMP (Fig. 3, panel i) and the dissociation of the MMP9⅐TIMP complex (Fig. 3, panel ii) against the experimental data. The kinetic parameters that we determined are tabulated in Table 2.
Schemes II and III were motivated by characterization of the 62-kDa MMP2 interaction with TIMP2 as shown by Hutton et al. (25). A similar mechanism was also previously considered between MMP1 and TIMP1 wherein the intermediate complex represented TIMP1 bound to the hemopexin domain of MMP1 (26). Scheme III is a variation of scheme II in that the intermediate is considered as fully active as the 82-kDa MMP9. Jointly, schemes II and III represent the extremes of the general binding ϩ isomerization model, which we consider as scheme III*. Scheme IV was suggested by O'Connell et al.

TABLE 1 Kinetic parameters of activation
The activation of pro-MMP9 by various activators was modeled using a Michaelis-Menten reaction scheme. The kinetic parameters are listed in the 1st column for the respective proenzyme and activators. The on-and off-rates of activator towards the various pro-MMP9⅐TIMP⅐MMPX complexes were assumed equal to that towards pro-MMP9. (f indicates fitted parameter; A indicates assumed parameter based on comparison to similar reactions.)

Kinetic parameters of activation at 37°C
(pro-species in parentheses) Ref. DECEMBER 28, 2007 • VOLUME 282 • NUMBER 52 (18) in noting two subpopulations of MMP9, with one capable of significantly stronger inhibition. We found that scheme I can accurately characterize the inhibition of 82-kDa MMP9 by TIMP2 (and of 62-kDa MMP2 by TIMP1; data not shown) ( Fig. 3A) but not by TIMP1 (Fig. 3B).

Biochemical Model of MMP9 Activation and Inhibition
The differences between MMP9 inhibition by TIMP1 and TIMP2 cannot be accounted for simply by altering the kinetic parameters used in a scheme I reaction mechanism. On the other hand, we found that schemes III* and IV seem capable of characterizing inhibition between both MMP9 and TIMP1 and The panels are as follows: MMP9 ϩ TIMP2 using scheme I (A) and MMP9 ϩ TIMP1 using schemes I (B), II (C), III (D), III* (E), and IV (F). Note: for MMP9 inhibition with TIMP2, we find that the experimental data in A, panel ii, displays an onset of activity that appears more linear than the scheme I produced fit. We found that schemes II and III can approximate this linearity; however, the improvement was deemed indicative of overfitting and was also found to be unimportant to later conclusions in Fig. 5.
MMP2 and TIMP2 (latter not shown). The results for TIMP1 can also be seen by comparing the normalized errors for each of the schemes in Table 2.
Scheme II (see Fig. 3C) does not significantly improve upon the results of scheme I and does not also seem capable of describing the MMP9/TIMP1 interaction. We searched the space of parameters for k 1 II between 1 ϫ 10 3 and 1 ϫ 10 8 M Ϫ1 s Ϫ1 , and for k 2 II , k 3 II , and k 4 II between 1 ϫ 10 Ϫ8 and 1 ϫ 10 2 s Ϫ1 ; the approximate error-minimizing fit is shown in Fig. 3C, using a k 1 II ϭ 2 ϫ 10 6 M Ϫ1 s Ϫ1 , k 2 II ϭ 4 ϫ 10 Ϫ3 s Ϫ1 , k 3 II ϭ 6 ϫ 10 Ϫ4 s Ϫ1 , and k 4 II ϭ 3 ϫ 10 Ϫ5 s Ϫ1 , with an overall K i of k 2 II /k 1 II ϫ k 4 II /(k 3 II ϩ k 4 II ) equal to 96 pM. This is expected because the initial binding provides inhibition, and the subsequent isomerization to another inhibited state can only increase the rate of the initial binding inhibition.
Scheme III and more so schemes III* and IV come closer to approximating the MMP9/TIMP1 interaction. In the case of scheme III (Fig. 3D), the model produced a biphasic approximation of the data; however, the proper active MMP9 concentration dependence on TIMP1 could not be reproduced. The effective K i ϭ k 4 III /k 3 III ϫ k 2 III /k 1 III of 1.4 pM is well in the range approximated by O'Connell et al. (18) at K i Ͻ 50 pM. Using scheme III* (Fig. 3E), we determined that setting the activity of the intermediate complex to ϳ35% that of free MMP9 provided optimal fitting for the isomerization type schemes. For scheme III*, the other parameters were determined at k 1 III* ϭ 5.6 ϫ 10 6 ) or 0.9 pM. Similar to scheme III*, scheme IV is also able to reproduce the biphasic behavior and the proper TIMP1 concentration dependence (Fig. 3F). However, this occurred only when the ratio of the slowly inhibited (low affinity MMP9 L ) to total MMP9 (MMP9 L ϩ MMP9 H ) species was set to ϳ0.31. The study of O'Connell et al. (18) determined that the rapidly inhibited MMP9 form had a k 1 IV ϭ 1.11 ϫ 10 7 M Ϫ1 s Ϫ1 . We estimated the parameters of this model to be k 1 IV ϭ 4.6 ϫ 10 6 M Ϫ1 s Ϫ1 and k 2 IV ϭ 6.3 ϫ 10 Ϫ6 s Ϫ1 for the rapidly inhibited form (K i ϭ 1.4 pM) and k 3 IV ϭ 2.9 ϫ 10 3 M Ϫ1 s Ϫ1 , k 4 IV ϭ 2.5 ϫ 10 Ϫ5 s Ϫ1 for the slowly inhibited form (K i ϭ 8.6 nM). We note that the rapidly inhibited form also was consistent with the strong inhibition strength seen by O'Connell et al. (18) and that the low affinity isoform in the range noted by Olson et al. (19) The previous estimates are the optimal sets of parameters produced by our grid-search algorithm. However, as with any nonlinear parameter estimation algorithm, we have no guarantee for finding the absolute "global" error minimizing set of rate parameters. Furthermore, although there may be only one best fit, there may be a family of nearly well fitting curves that cannot be distinguished both on the basis of experimental uncertainty and the inability of experimental data to confer values for parameters of a specified model. Thus, similar to our attempt in Fig. 2B, we attempted to specify ranges of parameters (more accurately, the region traversed by the error-minimizing manifold). We only did this for schemes III* and IV because they were the only models consistent with the experimental data. We found that the error-minimizing manifold for scheme IV traversed over the following ranges: The ratio of the low affinity isoform to total MMP9 was relatively defined at 0.3-0.35. For these parameters, the effective K i value for MMP9 H varied from 0.2 to 1.5 pM. Similarly for MMP9 L , the K i value ranged from 0.4 to 15 nM. The scheme III* error-minimizing manifold traverses a wider range of parameter combinations: We also tested the 62-kDa MMP2/TIMP2 interaction using the various inhibition schemes, again using data from Olson et al. (19). Our conclusions are similar to that found for the MMP9/TIMP1 interaction; schemes III, III*, and IV are much better at characterizing the data than schemes I and II. Previous experimental estimates by Hutton et al. (25) also characterize this reaction. They used the reactions of scheme II/III to estimate on and off rates for the initial association at 5.9 ϫ 10 6 M Ϫ1 s Ϫ1 and 6.3 s Ϫ1 , respectively, and forward and reverse isomerization rates at 33 and 2.0 ϫ 10 Ϫ8 s Ϫ1 , respectively (effective K i ϭ 0.6 fM) (25). These kinetic parameters were not able to simulate the behavior found in Olson et al. (19).
Although scheme IV is theoretically interesting, experimental work needs to be done to determine whether two isoforms of MMP9 do exist in biochemical preparations. For this reason, we used scheme III* to describe the MMP9/TIMP1 interaction in subsequent results. One possibility is that the two isoforms arise because of differences at the time of protein synthesis (e.g. glycosylation or folding). The difference is TIMP1-specific (because scheme IV is not needed to model the MMP9/TIMP2 interaction); it is probable that these differences are located in the MMP9 C-terminal domain. It is also tempting to speculate that post-secretion processes such as C-terminal truncation can also result in these TIMP1-specific differences. This is cor-

Mechanisms of MMP9 inhibition
We numerically fitted kinetic parameters for MMP9 inhibition by TIMP using schemes I-IV and data from Olson et al. (19). Reaction diagrams of the inhibition schemes are presented in Fig. 3. The chosen set of parameters was determined to be optimal by our grid-search algorithm. The normalized (Norm.) error measure used was 1 Ϫ R 2 . The ranges of parameter sets for schemes III* and IV are presented in the text.

Scheme
no.

Norm. error
Effective DECEMBER 28, 2007 • VOLUME 282 • NUMBER 52 JOURNAL OF BIOLOGICAL CHEMISTRY 37589 roborated by experimental studies indicating that MMP9 is subject to C-terminal processing by MMP3 (13), which was used in pro-MMP9 activation in the experimental setup of Olson et al. (19), and data indicating that C-terminal truncated forms of MMPs have reduced affinities to TIMPs (18,19,27).

Biochemical Model of MMP9 Activation and Inhibition
Based on scheme IV, we can also consider the influence of C-terminal processing on MMP9 inhibition by TIMP1. We found that varying the ratio of the two forms of MMP9 in scheme IV produced drastically different progress curves for inhibition by TIMP1 (not shown). This implies that if C-terminal processing is biochemically relevant, then the experimental determination of MMP9 inhibition by TIMP1 will also be very sensitive to the ratio between the 82-and 65-kDa MMP9 isoforms, i.e. to the precise concentrations of activators and duration over which MMP9 was activated.
Homolytic Processing Produces Second-order MMP9 Decay Kinetics-Thus far we have not included deactivation (loss of MMP9 activity because of conversion to ϳ50 -60 kDa by processing at the MMP9 active site) or C-terminal cleavage (conversion to 65-kDa MMP9) into the model as there have been conflicting results in the literature on whether these processes occur and how. Whether these processes are physiologically relevant is another question because any existing TIMPs will significantly delay these processes (15). We analyzed the existing data that support deactivation. We did not examine C-terminal processing because of a lack of quantifiable data.
The existing experimental observations show conflicting results in whether the C-terminally truncated and deactivated isoforms are present. For example, Okada et al. (13) observe a reduction in MMP9 activity leading to ϳ25% activity by 54 h upon activation by MMP3. SDS-PAGE reveals conversion to the 65-kDa isoform but no subsequent conversion to the ϳ50 -60-kDa isoform. Similar trends were observed with trypsin as an activator (see Fig. 4, A and B). This is in contrast to the results of Nakamura et al. (11) and Shapiro et al. (15) where activation by MMP10 and MMP3 clearly leads to formation of both the ϳ65and 50 -60-kDa isoforms. In contrast, SDS-PAGE results in O'Connell et al. (18) seem to indicate direct conversion from the 82-kDa MMP9 to an ϳ40-kDa isoform after activation by MMP3 for 22 h.
We make several clarifications in our treatment of MMP9 deactivation. Because we have relied on normalization of experimental data in our estimation of kinetic parameters for MMP9 activation, accounting for deactivation would affect our determination of these parameters (e.g. the maximum activity in the presence does not represent 100% activation). However, we found that the changes to the kinetics of MMP9 activation were not significant enough to be resolved by the experimental data. Additionally, deactivation should not necessarily be limited to MMP9. Mechanisms may exist that also allow the deactivation of the activating species, TIMPs, experimental substrates used for activity quantification, or even pro-MMP9. We attempt to consider these later. Finally, for simplicity we assume that deactivation occurs via a single step irreversible bimolecular reaction (instead of a Michaelis-Menten scheme as we have done for activation).
Based on the literature, we constructed several possible schemes of deactivation (Fig. 4). These include a heterolytic scheme in which the activator of pro-MMP9 catalyzes deactivation of MMP9, a bimolecular homolytic scheme in which an active MMP9 deactivates another active MMP9, and a unimolecular reaction scheme which can be used to represent the autolysis of a single MMP9 molecule or the loss of MMP9 activity because of unfavorable environment conditions.
To derive parameters for the deactivation kinetics, we adopted experimental data by Okada et al. (13) and Sang et al. (24) (Fig. 4). We only present the former, in whose study pro-MMP9 reacts with two different concentrations of trypsin (420 and 42 nM) (Fig. 4A) and MMP3 (167 and 33 nM) (Fig. 4B). This study found near equal rates for degradation at both concentrations of initial activator, despite the activator used.
Our simulations (Fig. 4) suggest that the bimolecular homolytic mode of deactivation is the only scheme consistent with the selected experimental data. This is in conflict with observations that the rate of deactivation is activator-specific (13,24). We did find dependence of the rate constant of homolytic degradation, k hom , on the activator used (e.g. 1300 and 450 M Ϫ1 s Ϫ1 for activation by trypsin and MMP3, respectively); however, within each experiment, the homolytic mode of deactivation produces the most consistent results. The homolytic rates of deactivation are approximately 1 order of magnitude less than the k act /K m values for MMP9 activation. The autolytic scheme is also able to display reasonable fit for k auto , the autolytic rate constant, between 1.0 ϫ 10 Ϫ5 and 1.5 ϫ 10 Ϫ5 s Ϫ1 in the presence of MMP3 and around 2 ϫ 10 Ϫ5 s Ϫ1 in the presence of trypsin. In contrast, the heterolytic scheme is completely inconsistent with data for both MMP3 and trypsin, and the optimal rate constant for heterolytic deactivation was found to be 0 M Ϫ1 s Ϫ1 . The three schemes reflect differences of the time course of MMP9 decay in both the rate of decay of activity and the dependence on activator concentration. The autolytic and heterolytic schemes both exhibit very rapid deactivation (after sufficient MMP9 is activated) that can be characterized by a firstorder "exponential" decay (see Fig. 4B). On the other hand, the experimental data display a slower rate of deactivation, which can be characterized by the second-order "hyperbolic" time dependence of homolytic deactivation. In terms of dependence on the activator concentration, the experimental data show both a relative insensitivity to changes in the activator concentration and positive correlation between the activator concentration and peak MMP9 activity. The autolytic and homolytic deactivation schemes reflect this behavior. However, the opposite is true for the simulated heterolytic scheme. From Fig. 4, A and B, we see that the heterolytic scheme results in normalized activity greater than 100% at the lower concentration of activator. This surprising behavior occurs because of the fact that decreasing the activator concentration decreases the rate of MMP9 activation, which further adds to the decreased rate of MMP9 deactivation. This results in greater MMP9 accumulation at lower activator concentrations.
Other mechanisms may also be at work. The presence of nonspecific activity would result in final levels of activity greater than zero. This would give the false appearance of a slower approach to the final equilibrium, suggesting that the observed deactivation may indeed be autolytic and a result of poor storage conditions. However, this does not appear to be the case as the activity of MMP3 was also measured in Ref. 13 and was shown to be negligible. Next, although we showed that the overall decay was second order, we assumed this was entirely the result of a second-order decay in MMP9. If MMP9 was reacted in sub-optimal storage conditions, other alternatives can be present yielding possibly overall second-order decay via the following: 1) both the quantification substrate and MMP9 decaying with first-order kinetics; 2) both MMP3 (or trypsin) and MMP9 decaying with first-order kinetics in addition to heterolytic deactivation of MMP9 by MMP3; and 3) both pro-MMP9 and MMP9 autolytically decaying with first-order kinetics. We tested these possibilities, but the first and third are only marginally better than the autolytic reaction and still result in exponential decay at the longer times of this experiment. An interesting possibility is the second case. Deactivation of both MMP3 and MMP9 with a strong autolytic deactivation for MMP3 and a weak autolytic and strong heterolytic deactivation for MMP9 does produce second-order decay kinetics in MMP9; however, it still displays the same concentration dependence as the original heterolytic scheme. Thus, the data we used are only consistent with a homolytic bimolecular decay mechanism for MMP9.
TIMPs Delay the Onset of MMP9 Activation; Importance of the MMP:TIMP Ratio-Previous biochemical experiments indicate that pro-MMP9 activation is severely limited when TIMP concentrations exceed that of a matrix metalloproteinase pro-MMP9 activator. To determine how the competitive inhibition by TIMPs influences the onset of MMP9 activation, we sought to characterize the dynamics of the MMP3/TIMP2/ pro-MMP9 system. We did not take into account deactivation of MMP9.
We simulated activation (Fig. 5, A and B) of 300 nM pro-MMP9 by 50 nM MMP3 in the presence of 0, 20, 40, 50, 60, 75, and 100 nM TIMP2. As expected, the time required for complete MMP9 activation is increased as TIMP2 concentration is increased. When TIMP2 reacts stoichiometrically with MMP3, both activity (Fig. 5A) and activation (Fig. 5B) profiles develop an initial concavity, which increases in magnitude as the inhibitor concentration is increased. At excess inhibitor, activation profiles reveal a steady but greatly reduced rate of activation, at which time MMP9 activity is nearly completely inhibited. This trend is easily visible at 100 nM TIMP2, where activity is delayed by 18 h. These distinct behaviors were present because our concentrations were well above the K i of the inhibition reactions. In our Supplemental Material, we tested a larger range of MMP3 and TIMP2 concentrations to determine the effects of TIMP on the delay in MMP9 activation, activity, and deactivation.
We can dissect the mechanism of the lag in MMP9 activation and activity by computing the concentrations of the major species (Fig. 5C) and the velocities of the reactions forming MMP9 (Fig. 5D) for the case of 100 nM TIMP2. After an initial relaxation between MMP3, TIMP2, and pro-MMP9 (Fig. 5D before 1 h), the rate of production of MMP9 by MMP3 is almost equally balanced by the rate of inhibition of MMP9 by TIMP2 for up to 18 h. This is because any MMP9 that is formed is rapidly inhibited by the high concentration of TIMP2. Thus, the lag in activity is because of the low levels of free MMP3 and the high levels of free TIMP2. In this period, the overall [MMP]: [TIMP] (measuring all activated and complexed species) is between 1 ⁄ 2 and 1. At 18 -22 h, the depletion of TIMP2 by MMP9 (resulting in the dissociation of the MMP3⅐TIMP2 complex) frees MMP3 and results in an acceleration of MMP9 activation. Thus, MMP9 activity sharply rises in this period as the total MMP overtakes total TIMPs, i.e. when the overall [MMP]: [TIMP] surpasses 1.
Experimentally, many studies show lack of MMP9 activation for very long durations in conditions where TIMP1 exceeds MMP3 (15,28). In Hahn-Dantona et al. (28) , an excess of 1.5 nM TIMP1 over 2.5 nM MMP3 leads to no significant production of the 82-kDa MMP9, until up to 24 h after mixing. However, our results show that MMP9 should be slowly but steadily activated (but inactive because of inhibition by TIMP) because a small fraction of MMP3 should exist in its free state. The fact that activation is not readily observed may be due in part to a tighter than expected MMP3⅐TIMP1 complex or to some sufficient level of deactivation.
Pro-MMP9 Activation by MMP3 in the Presence of TIMP1 Possibly Leads to Overshoot in MMP9 Activity-Next, we attempted to extend the model to encompass pro-MMP9 activation by MMP3 and trypsin in the presence of TIMP1. In this system, pro-MMP9 can bind TIMP1 leading to the pro-MMP9⅐TIMP1 and pro-MMP9⅐TIMP1⅐MMP3 intermediates. These complexes can remain intact through activation leading to the active MMP9⅐TIMP1* and MMP9⅐TIMP1⅐MMP3 complexes (29). We consider that MMP9⅐TIMP1*, which is likely a different intermediate than that formed in the scheme III MMP9/TIMP1 reaction, can undergo reactions similar to the latter. Uncertainties in our model include the rates of activation (k on , k off , and k act ) of the proenzyme intermediates and the level of activity and the modes of dissociation of the activated complexes. As possibilities, we considered that the activated complexes could dissociate at either interaction site (see legend to Fig. 6, D and E) or dissociate instantaneously following activation (see legend to Fig. 6F). We have attempted to shed light on these unknowns by attempting to simulate experimental data from three separate studies (Fig. 6, A-C) (28 -30). Because of the increased complexity of the system, we make several simplifications.
First, we assume that activation of the pro-MMP9⅐TIMP1 and pro-MMP9⅐TIMP1⅐MMP3 complexes by MMP3 occurs with the same kinetic on and off rates. We lump any differences in the net activation rate into k act . Specifically, we expect that the values of k act be smaller for the complexed intermediates than for free pro-MMP9 because of steric hindrance. We also assume that the reduction in k act is the same for both complexes.
Second, pro-MMP9⅐TIMP1 inhibits MMP3 to form pro-MMP9⅐TIMP1⅐MMP3. Ogata et al. (30) have suggested that there is increased dissociation between pro-MMP9 and TIMP1 when this happens. In contrast, Kolkenbrock et al. (29) have shown remarkably stable MMP9⅐TIMP1⅐MMP3 complexes. Based on our analysis in the Supplemental Material (see Determining Reaction Cycles), we expect either a strongly bound or weakly bound pro-MMP9⅐TIMP1⅐MMP3 complex. We consider the former case because the latter would resemble the MMP3/TIMP2/pro-MMP9 system. Thus, we would expect activation to occur mainly through the pro-MMP9⅐TIMP1⅐MMP3 complex.
Finally, we can hypothesize that pro-MMP9⅐TIMP1 can also inhibit active MMP9 just as it can MMP3. However, experimental results have demonstrated that MMP9 is weakly inhibited by pro-MMP9⅐TIMP1 despite its strong capacity to be inhibited by TIMP1 alone (29). Assuming weak binding kinetics between MMP9 and pro-MMP9⅐TIMP1, we have verified that this reaction is negligible (data not shown); thus, we do not consider it further.
The experimental results (28 -30) and our simulations are shown in Fig. 6, A-C. Fig. 6A describes MMP9 activity following activation of pro-MMP9 by increasing levels of MMP3 in either the presence or absence of 4 nM TIMP1 (28). Similar to the experimental data, we found that activation of MMP9 only occurs for MMP3 in excess of TIMP1. In Fig. 6B, the pro-MMP9⅐TIMP1 complex was activated by trypsin for varying amounts of MMP3 (29). This study also showed that without trypsin, MMP3 concentration must first exceed that of the pro-MMP9⅐TIMP1 complex for MMP9 activity to ensue. Another interesting result from this study was the partial activity of the trypsin-activated pro-MMP9⅐TIMP1 complex (see Fig. 6B, top curve at added MMP3 ϭ 0 ng). The final experiment we considered (Fig. 6C) was adopted from Ogata et al. (30) and monitored the activity of MMP9 at 1 and 4 h following reaction of pro-MMP9⅐TIMP1 by varying concentrations of MMP3. All three experiments showed that MMP3 needed to exceed the concentration of TIMP1 for emergence of MMP9 activity. The reaction mechanisms and kinetic parameters we found most suitable are presented in Fig. 1 and Tables 2 and 3.
We find that the simulated results provide a reasonable approximation of the experiments (Fig. 6, A-C). However, a significant discrepancy is visible in Fig. 6C However, our calculations produce significantly greater activity at these points, which we found was a result of dissociated MMP3 activating pro-MMP9. This indicated that the experimental data require a stronger inhibition of MMP3 by TIMP1 and pro-MMP9⅐TIMP1. However, to numerically achieve as strong an MMP3 inhibition as possible, we have assumed that pro-MMP9⅐TIMP1 inhibits MMP3 with the same kinetics as free TIMP1.
The different modes of dissociation of the activated pro-MMP9⅐TIMP1 and pro-MMP9⅐TIMP1⅐MMP3 complexes that we considered (see Fig. 6, D-F, legend) were all capable of reasonably reproducing the observed data, albeit for slightly different k act values. Despite the similar appearance against the experimental data, analysis of the time courses of the different mechanisms revealed significantly different temporal phenomena (Fig. 6, D-F). When we considered that MMP9⅐TIMP1* and MMP9⅐TIMP1⅐MMP3 intermediates were active, we found that MMP9 activity initially overshot its equilibrium level to return to it at later times (Fig. 6, D and E). This behavior has not been previously characterized for the MMP/TIMP systems. As presented in the Supplemental Material, the overshoot occurs because the rate of MMP9 activation exceeds the rate of MMP9 inhibition. This also implies a role for pro-MMP9/TIMP1 complexation in the persistence of overshoot. On the other hand, if we assume that either the activity of the activated complexes is limited or that the complexes dissociate rapidly to become instantly inhibited in MMP9, then little or no overshoot occurs (Fig. 6F). We found that the presence of overshoot in this system was very sensitive to rates of MMP3 inhibition, MMP9 inhibition, and MMP9 activation and to the structure of the reactions. We did not find overshoot behavior in the MMP3/ TIMP2/pro-MMP9 system unless the kinetic parameters of the system were significantly altered (e.g. greater than a 10-fold reduction in rate of inhibition; data not shown).
In the MMP3/TIMP1/pro-MMP9 system, our simulations only indicate that overshoot is a possibility. The experimental evidence that supports the notion of overshoot in the activation of pro-MMP9⅐TIMP1 by MMP3 or trypsin is limited (21,28). In contrast, activation of pro-MMP9⅐TIMP1 by trypsin by Ogata et al. (30) reveals no noticeable overshoot. In addition, what is sometimes observed as activation followed by deactivation may in fact be activation followed by inhibition.
An implication of this result is that the experimental data for MMP9 activity in Fig. 6C does not necessarily have to be monotonically increasing in time, as shown in Fig. 6, D and E. To achieve monotonicity in MMP9 activity between 1 and 4 h, we had to delay the period of overshoot by lowering the effective

Parameters used in MMP interactions with TIMP1 and TIMP2
We tabulated the kinetic parameters used in our simulations for the pro-MMP9/ TIMP1, MMP9/MMP3⅐TIMP1, MMP3/TIMP1, and MMP3/TIMP2 interactions. Parameters were taken from the experimental literature. k off between MMP3 and TIMP2 was assumed equal to that between MMP3 and TIMP1. k act of the pro-MMP9⅐TIMP1 and pro-MMP9⅐TIMP1⅐MMP3 complexes to less than 40% of the k act of pro-MMP9 alone. However, as this result is possibly a numerical artifact of our assumptions, further biochemical studies would be required to determine the significance of overshoot in pro-MMP9⅐TIMP1 activation.
Our simulation of the experimental procedure by Kolkenbrock et al. (29) (Fig. 6B) also suggests a different interpretation of its experimental results. In Fig. 6B, activation of pro-MMP9⅐TIMP1 using trypsin and zero MMP3 results in some nonzero level of activity. However, the addition MMP3 (up until an initial concentration equal to that of pro-MMP9⅐TIMP1) subsequently results in increased MMP9 activity. The original conclusion that activation of pro-MMP9⅐TIMP1 (at zero added MMP3) by trypsin yields a complex with less activity than activation of pro-MMP9⅐ TIMP1⅐MMP3 (for added MMP3) does not factor in the ability of MMP3 to sequester TIMP1, effectively freeing MMP9 from inhibition. Thus numerically both MMP9⅐TIMP1* and MMP9⅐TIMP1⅐MMP3 can have equal activity and still provide simulation of experimental data.
Conclusion-Construction of a biochemically detailed model of MMP9 allowed us to revise and extend conclusions drawn from the previous literature. We were also able to provide kinetic parameters for MMP9 activation, inhibition, and deactivation. We summarize our primary findings and predictions in Table 4.
We should note that in our parameter estimation, we often found a range of acceptable sets of kinetic parameters that allowed simulation of experimental data. Generally, this results from fitting a complex model with experimental data that cannot directly probe certain areas of the "kinetic parameter space" of the model. However, our results were often insensitive to the exact parameters we used in our models. For example, if two different sets of kinetic parameters for MMP9 inhibition produced similar fits against progress curves of inhibition, then the results remained similar when in the presence of pro-MMP9 activators. Additional experiments would need to be undertaken to accurately confer values for the degenerate parame-ters. Another problem typical in estimation of parameters is that of over-fitting complex models. However, as all of our schemes are biochemically reasonable and as we generally picked the simplest scheme capable of providing "visual" accuracy against experimental data, we do not expect severe overfitting.
Regarding MMP inhibition, we found two mechanisms equally capable of describing biphasic inhibition. Although Hutton et al. (25) described a scheme similar to scheme III*, we also found that scheme IV is capable of producing progress curves exhibiting biphasic inhibition. Overall, our simulations confirm previous observations regarding the origins of this inhibition (18,25). An interesting correlation is that both MMP9/TIMP1 and MMP2/TIMP2 form proenzyme inhibitor complexes. It would be interesting to determine whether the temporal aspects of biphasic inhibition also are present between MMP9/TIMP3, MMP2/TIMP3, and MMP2/TIMP4, which are also capable of proenzyme inhibitor complexes. This may help to further clarify the role of the C-terminal domain and answer which of schemes III* and IV is involved in the observed dynamics of inhibition. Evidence for scheme IV would come from simultaneous direct observation of two or more isoforms of MMP9 and a dependence of the temporal dynamics of inhibition (at a fixed active MMP9 concentration) on the protocol used for MMP9 activation.
We should also note that despite using the progress curve data of Olson et al. (19), our determination of MMP9 inhibition by TIMP1 produced K i values that were at least 3 orders of magnitude smaller, consistent with the pM range of K i found previously (18,20,33). This may be a result of different methodology. Whereas Olson et al. (19) (and other experimental studies) have used estimates from kinetic theory derived for slow binding inhibitors (34), we have numerically integrated the differential equations resulting from our various reaction mechanisms. Both methods depend on the assumptions used between enzyme, inhibitor, and substrate binding. In our method, we have no guarantee of finding the globally errorminimizing set of parameters. However, we were able to come to within a factor of 1.3 of the error-minimizing manifold; thus Pro-MMP9 activation in conditions where TIMP exceeds an MMP activator will display a lag in MMP9 activity followed by a relatively sharp rise in the rate of activation No 4 The intermediates of pro-MMP9⅐TIMP1 and pro-MMP9⅐TIMP1⅐MMP3 complex activation are active, but the relative activities of the two to TIMP-free MMP9 remain unknown No 5 Pro-MMP9⅐TIMP1 complexes are activated at a reduced rate compared with free MMP9 No 6 Pro-MMP9 TIMP1 complexation may provide two competing effects on MMP9 activity. It may result in reduced rates of activation of complexes, but it may also lead to increased persistence of the active forms Yes 7 The MMP3/TIMP2/pro-MMP9 system does not exhibit an overshoot in MMP9 activity, whereas the MMP3/TIMP1/pro-MMP9 system does Yes we do not expect our methodology to account for the large difference in magnitude of K i . At present, it is not clear what other factors may have led to the different results. Next, we studied MMP9 deactivation, a process for which little is known. Most of the kinetic studies we have encountered showed some degree of deactivation. Whether this is an artifact of experimental protocol or a feature of the MMP and activator system can only be revealed by more accurate quantification and reporting. Our simulations support homolytic deactivation. This result is contrary to the typical assumption of heterolytic deactivation when deactivation does occur. However, the experimental data we studied, Okada et al. (13) and Sang et al. (24), clearly show that this process cannot be purely activator-driven, on the grounds that MMP9 deactivation seemed largely insensitive to the activator concentration. Because of significant differences in deactivation between different studies, more detailed experimental work in this field is definitely needed.
Finally, our attempt to characterize MMP9 activation by MMP3 in the presence of MMP inhibitors revealed a novel temporal behavior where the activity of MMP9 increases because of activation but then decreases because of a slower inhibition response. Despite TIMP1 being the stronger inhibitor of MMP9 than TIMP2, we found this overshoot behavior prevalent in our numerical simulations (and sensitivity analyses) of the MMP3/TIMP1/pro-MMP9 system, but not in those of the MMP3/TIMP2/pro-MMP9 system. This overshoot is enabled when the process of inhibition occurs in series with the process of activation (as opposed to in parallel). In our case, overshoot was mainly a result of the accumulation of the MMP9⅐TIMP1⅐MMP3 intermediate, which was assumed to be incapable of inhibition. However, this assumption might be accurate as quaternary complexes are not observed in the literature. On the other hand, this complex has also been shown to be capable of substrate proteolysis, which could have also resulted from partial dissociation of the MMP9. If we naturally assume the activity of the complex is diminished compared with free MMP9, the overshoot in MMP9 activity does indeed decrease and is delayed. We found similar behavior when we assumed decreased rates of formation or increased rates of dissociation of the intermediates (or even a weakened MMP3/TIMP1 interaction allowing more TIMP1 to inhibit MMP9). Underlying these is the balance between the rates of MMP9 formation and inhibition.
The data we analyzed can support overshoot, but this is not necessary (Fig. 6, D-F). Direct experimental data showing overshoot may also exist (14); however, they need to be confirmed in the absence of MMP9 deactivation. Ultimately, this can lead to simplification of complex systems such as MMP3/TIMP1/pro-MMP9 into the "effective" simpler systems such as MMP3/TIMP2/pro-MMP9.
The cellular context in which the studied processes are relevant is still being explored. These conclusions will hopefully help guide future efforts to understand complex proteolytic systems, which we are pursuing in the context of angiogenesis. In this arena, MMP9 is hypothesized to play a variety of functions from proteolyzing type IV collagen in basement membranes surrounding activated endothelial cells to releasing vascular endothelial growth factor sequestered on glycosaminoglycan chains of proteoglycans. Using computational models, we can numerically study the effect of MMP9 and other MMPs and TIMPs on growth factor-mediated angiogenesis. Ultimately, MMP9 and TIMPs may have other yet unknown functions. However, the importance of computational models lies in their ability to verify if current biochemical understanding is physically realizable, especially in dynamic and spatially complex cellular environments.