Cycle Network Model of Prostaglandin H Synthase-1

The kinetic model of Prostaglandin H Synthase-1 (PGHS-1) was developed to investigate its complex network kinetics and non-steroidal anti-inflammatory drugs (NSAIDs) efficacy in different in vitro and in vivo conditions. To correctly describe the complex mechanism of PGHS-1 catalysis, we developed a microscopic approach to modelling of intricate network dynamics of 35 intraenzyme reactions among 24 intermediate states of the enzyme. The developed model quantitatively describes interconnection between cyclooxygenase and peroxidase enzyme activities; substrate (arachidonic acid, AA) and reducing cosubstrate competitive consumption; enzyme self-inactivation; autocatalytic role of AA; enzyme activation threshold; and synthesis of intermediate prostaglandin G2 (PGG2) and final prostaglandin H2 (PGH2) products under wide experimental conditions. In the paper, we provide a detailed description of the enzyme catalytic cycle, model calibration based on a series of in vitro kinetic data, and model validation using experimental data on the regulatory properties of PGHS-1. The validated model of PGHS-1 with a unified set of kinetic parameters is applicable for in silico screening and prediction of the inhibition effects of NSAIDs and their combination on the balance of pro-thrombotic (thromboxane) and anti-thrombotic (prostacyclin) prostaglandin biosynthesis in platelets and endothelial cells expressing PGHS-1.


Introduction
Prostaglandin H synthase-1 (PGHS-1), also known as Cyclooxygenase-1 (COX-1), catalyses the first stage in the conversion of arachidonic acid (AA) into prostaglandin H 2 (PGH 2 ), being a precursor of physiologically active prostanoids [1]. PGHS-1 and the other isoform of PGHS (PGHS-2) are the main targets of non-steroidal anti-inflammatory drugs (NSAIDs) in the treatment of pain, fever, inflammation, and cardiovascular diseases [2]. Despite the fact that this enzyme has been under thorough investigation for several decades, some details of its catalytic mechanism are still not completely understood [3]. It is generally recognized that PGHS-1 has two distinct catalytic activities: a cyclooxygenase (COX) activity converting AA to an intermediate product-prostaglandin G 2 (PGG 2 )-and a peroxidase (POX) activity which further converts PGG 2 into a final product-PGH 2 . PGHS-1 catalysis represents a complex network of many oxidized/reduced reactions connecting multiple enzyme intermediates and leading to redox transformation of key catalytic components such as the arachidonic acid (AA) binding site (COX site) and the heme-prosthetic group (POX site). Another not completely understood phenomenon in PGHS-1 catalysis is the so-called self-inactivation process, including gradual decay of both COX and POX activities during PGHS-1 catalysis [4]. Commonly, two different catalytic Figure 1. Scheme of the branched chain mechanism of Prostaglandin H synthase-1(PGHS-1) catalysis [5]. Detailed description of all the enzyme states is given in Table 1.
The most detailed mathematical models of PGHS-1 were developed in a series of papers [5,[8][9][10]. These kinetics models were based on the assumptions of the BC mechanism of PGHS-1 function, and gradually evolved starting from eight reactions in an initial model [8] to 22 reactions in the latest version [10]. An increase in the complexity of the models reflected a growing insight into the mechanism of catalysis and molecular structure of PGHS-1. An extension of the models included an introduction of additional enzyme intermediates and reactions into these proposed catalytic schemes and more detailed consideration of enzyme self-inactivation. Each step in a gradual refinement of the model and an increase in its complexity were supported by the new experimental data as they became available. The developed models allowed satisfactorily describing key features of PGHS-1 catalysis, Figure 1. Scheme of the branched chain mechanism of Prostaglandin H synthase-1(PGHS-1) catalysis [5]. Detailed description of all the enzyme states is given in Table 1.
The most detailed mathematical models of PGHS-1 were developed in a series of papers [5,[8][9][10]. These kinetics models were based on the assumptions of the BC mechanism of PGHS-1 function, and gradually evolved starting from eight reactions in an initial model [8] to 22 reactions in the latest version [10]. An increase in the complexity of the models reflected a growing insight into the mechanism of catalysis and molecular structure of PGHS-1. An extension of the models included an introduction of additional enzyme intermediates and reactions into these proposed catalytic schemes and more detailed consideration of enzyme self-inactivation. Each step in a gradual refinement of the model and an increase in its complexity were supported by the new experimental data as they Pharmaceuticals 2020, 13, 265 3 of 21 became available. The developed models allowed satisfactorily describing key features of PGHS-1 catalysis, i.e., interaction of POX and COX activities, self-inactivation and autocatalytic effect, as well as estimating some kinetic constants of the PGHS-1 catalytic reactions [11].
The abundance of accumulated experimental data on the PGHS-1 kinetics stimulated the development of various kinetic models of PGHS-1. Several PGHS-1 models have been developed with a view to describe individual datasets and explain certain peculiarities of PGHS-1 functioning. For example, the kinetic models of PGHS-1 developed by Vrzheshch et al. [12] allowed for qualitatively describing sets of experimental data on AA, O 2 and RC consumption under specific in vitro conditions. The kinetic models of two isoforms, PGHS-1 and PGHS-2, were developed to explain the differences in their kinetics and activation thresholds [13].
None of the existing mathematical models of PGHS-1 can provide quantitative descriptions of all key regulatory effects within a unified set of kinetic parameters. In other words, different models and sets of parameters should be used to accurately reproduce a wide spectrum of experimental data on PGHS-1 kinetics and predict its function in vivo in health and pathological conditions. It complicated significantly the use of the existing models of PGHS-1 for the in silico study of the inhibition effect of NSAIDs in different in vitro and in vivo conditions and for their application in the development of quantitative systems pharmacology (QSP) models, where prostaglandins pathway is one of the drug targets.
In spite of great progress in understanding PGHS-1 functioning, the problem of quantitative description of the main experimental facts on its reaction mechanism as well as the reproduction of experimental data obtained in different in vitro and in vivo conditions in the framework of a unified kinetic model still remains. The new experimental data adds new knowledge regarding the enzyme structures and functions which can play a role in prostaglandin biosynthesis [14].
Basing on the ideas of the BC mechanism [5] and their realisation in the kinetic models [8][9][10], in this paper, we developed the microscopic kinetic approach to model the complex network dynamics of intra-enzymatic reactions among multiple catalytically active intermediate states of the enzyme and reproduce the detailed kinetic mechanisms of the PGHS-1 function. Note that this model has been previously applied to the in silico description of NSAID action alone and in combinations and showed satisfactory predictive results [15,16]. In particular, the model was applied to the investigation of dose dependence of PGHS-1 for aspirin and selective NSAID, celecoxib, on concentrations of arachidonic acid, reducing cosubstrate and hydrogen peroxide, which were reported to differ dramatically in purified and cellular assays, as well as in in vivo conditions in different cells and organs [17]. The model was also applied to the analysis of decreasing aspirin efficacy at its administration in combination with selective NSAIDs [15,16]. However, many aspects of the model and features of the interaction of NSAIDs and PGHS-1 have remained beyond the scope of these publications. In this paper, we give a full description of the model assumptions, method of its parametrisation, validation, and model application to the analysis of experimental data on regulatory properties of PGHS-1. One of the main aims of this work is the construction of reliable and validated kinetic model of PGHS-1, which would incorporate current theoretical conceptions on the PGHS-1 function and would be capable to quantitatively reproduce much of the experimental data on prostaglandin biosynthesis in the framework of a unified set of kinetic parameters.

Development of the Cycle Network Model of PGHS-1
The derivation of the model of PGHS-1 catalysis represents a special task, as both POX and COX activities are known to be located in different enzyme domains within the same protein complex, and these activities can proceed independently of each other [18]. Another important feature of PGHS-1 catalysis is that it involves intramolecular electron transfer and formation of multiple redox and binding states of the catalytic domain. The start point of our modelling was the BC mechanism [5,9] Pharmaceuticals 2020, 13, 265 4 of 21 schematically represented in Figure 1. The full catalytic cycle of PGHS-1 includes two distinct catalytic activities: a cyclooxygenase (COX) activity converting arachidonic acid (AA) to an intermediate product (prostaglandin G 2 (PGG 2 )); and a peroxidase (POX) activity, which further converts PGG 2 into final product PGH 2 . Following the approach to the PGHS-1 modelling [9], we assumed that PGHS-1 catalysis could be presented as a network of microscopic reactions, each of them resulting in a corresponding change in the enzyme catalytic state. Under a catalytic state, we understand the state of the PGHS-1 catalytic domain, which in its turn is defined by a unique combination of the states of its key elementary catalytic constituents, located in the COX and POX active sites. In our model, we considered the following key catalytically active components: fatty acid binding site (COX site), heme prosthetic group (POX site), and Tyr385 residue in the COX site.
We defined the COX site microstates through the presence/absence of AA in the hydrophobic channel of the COX region. The POX state was determined by the redox states of heme iron and protoporphyrin group. Tyr385 residue was assumed to exist in either a ground or radical state. In addition, we considered states of fully inactive enzyme (FIE) [19]. The detailed description of all-possible states for each catalytic component, considered in our model, is given in Table 1. Figure 2 shows a kinetic scheme for the developed catalytic cycle of PGHS-1. The scheme is mainly based on the assumptions of the BC mechanism [5] and its further extension [9]. Additionally, we also included in the model some reactions corresponding to an early version of the TC mechanism of PGHS-1 catalysis [6]. Thus, our scheme represents a generalisation of two previously discussed catalytic mechanisms, and the results obtained in the framework of the model allowed us to draw the conclusions on the role of each mechanism in PGHS-1 catalysis. converts PGG2 into final product PGH2. Following the approach to the PGHS-1 modelling [9], we assumed that PGHS-1 catalysis could be presented as a network of microscopic reactions, each of them resulting in a corresponding change in the enzyme catalytic state. Under a catalytic state, we understand the state of the PGHS-1 catalytic domain, which in its turn is defined by a unique combination of the states of its key elementary catalytic constituents, located in the COX and POX active sites. In our model, we considered the following key catalytically active components: fatty acid binding site (COX site), heme prosthetic group (POX site), and Tyr385 residue in the COX site. We defined the COX site microstates through the presence/absence of AA in the hydrophobic channel of the COX region. The POX state was determined by the redox states of heme iron and protoporphyrin group. Tyr385 residue was assumed to exist in either a ground or radical state. In addition, we considered states of fully inactive enzyme (FIE) [19]. The detailed description of allpossible states for each catalytic component, considered in our model, is given in Table 1. Figure 2 shows a kinetic scheme for the developed catalytic cycle of PGHS-1. The scheme is mainly based on the assumptions of the BC mechanism [5] and its further extension [9]. Additionally, we also included in the model some reactions corresponding to an early version of the TC mechanism of PGHS-1 catalysis [6]. Thus, our scheme represents a generalisation of two previously discussed catalytic mechanisms, and the results obtained in the framework of the model allowed us to draw the conclusions on the role of each mechanism in PGHS-1 catalysis.  Table 1. Enzyme states marked by colour correspond to the branched chain mechanism [5].  Table 1. Enzyme states marked by colour correspond to the branched chain mechanism [5].  We have presented the full catalytic cycle of PGHS-1 as a network of eight interlinked cycles-five of them (POX i ) corresponding to the POX activity, and the other three (COX j ) to the COX activity of the enzyme. In general, the full catalytic cycle of PGHS-1 can be considered as a complex dynamic network of elementary reactions which describes the totality of possible transitions between multiple states of the PGHS-1 catalytic domain.
Two interconnected cycles-COX 1 (reactions 1-4) and POX 1 (reactions 9-13)-correspond to the BC mechanism. According to this mechanism, PGHS-1 catalysis is initiated in the POX site, and its first step is the reaction of the resting enzyme state [Fe(III), PP] (E 1 ) with PGG 2 producing PGH 2 and leaving the heme prosthetic group in the state [Fe(IV), PP*+] (E 2 ). The enzyme state containing heme with two oxidizing equivalents-one on the iron and another one on the porphyrin radical cation-is referred to as Intermediate 1 [5]. In our model, we described this stage of POX catalysis by bimolecular irreversible reaction (reaction 11) of the resting state with PGG 2 , resulting in production of PGH 2 and leaving the enzyme in the E 2 state with heme in [Fe(IV), PP*+] state. The E 2 state can be converted back to the resting state E 1 by 2 consecutive reactions with RC by followed release of oxidized cosubstrate, OC (reactions 12 and 10). Alternatively, the E 2 state can undergo one-electron internal reduction via intramolecular electron transfer from Tyr385 residue in the COX site to protoporphyrin (reaction 13) resulting in the formation of the E 5 state [Fe(IV), PP, Tyr*] containing Tyr385 radical (Tyr385*) and a ferryl heme. This state is referred to as Intermediate II [5]. The E 5 state can also complete POX cycle and return to the resting state E 1 via reactions 9 and 10 with RC through the formation of [Fe(IV), PP] referred to as Intermediate II [5]. According to the BC mechanism, the presence of Tyr385 radical is ultimately required for AA oxidation in the COX site. Therefore, the intermediate state E 5 can potentially serve as a starting point for the further cascade of COX reactions.
The majority of previously known models included COX activity as a single Michaelis-Menten type reaction without considering its mechanism in any detail [10,13]. In contrast, in our model we describe COX activity by four consecutive reactions (reactions 1-4 in COX 1 cycle) corresponding to the transitions of the enzyme through three intermediate enzyme-substrate complexes (E 9 , E 13 , and E 8 ). Note that using Michaelis-Menten approximation to describe the individual reactions involved in a reaction network is not strongly valid.
Reaction 1 describes reversible binding of AA to the enzyme state E 5 , resulting in the formation of the state [Fe(IV), PP, Tyr*, AA] (E 9 ), containing AA in the COX site. Reaction 2 corresponds to the abstraction of a hydrogen atom from AA by Tyr385 radical and formation of the E 13 state containing AA radical in the COX site, [Fe(IV), PP, Tyr, AA*]. Further, AA * reacts with two oxygen molecules, and the intermediate product PGG 2 radical (PGG 2 *) is produced (reaction 3). We assumed that binding of the first and second oxygen molecules is a rapid process, during which the redox state of the heme remains unchanged [8]. The final step of the COX 1 cycle in our model is irreversible reaction 4, which describes the process of regeneration of Intermediate II with Tyr385 radical (E 5 ), accompanied by the release of the intermediate product PGG 2 from the COX site. Taking into account that POX activity of PGHS-1 is observed to proceed independently from COX catalysis [18], we considered the possibility of several additional POX and COX catalytic cycles, originating from intermediate states of the POX 1 and COX 1 cycles. For example, in the COX 1 cycle, the heme prosthetic group remains in the state [Fe(IV), PP] with one oxidising equivalent on iron, and therefore, each of the enzyme states E 5 , E 8 , E 9 , and E 13 can potentially participate in redox reactions with RC (reactions 19, 22, 28, and 29), producing catalytic states E 12 , E 15 , E 16 , and E 20 . The later states can be considered as intermediates of an additional COX 2 cycle. Note that similar COX cycle was included in the model developed by Kulmacz et al. [18]. The COX 2 cycle differs from the COX 1 cycle only by the state of the heme prosthetic group in the enzyme intermediates. Similarly, each intermediate of the POX 1 cycle has an empty COX binding site which allows for the binding of AA (reactions 30, 31, and 32) and production of the states E 6 , E 10 , and E 14 , respectively. These states can participate in the POX 2 cycle, which differs from the POX 1 cycle by the presence of AA in the COX binding site of its enzyme intermediates. The binding of AA to the resting state (E 1 ), containing Tyr385 in the ground state, was considered in the PGHS-1 model for the first time, as all previously known models assumed, that AA binds solely to the enzyme containing Tyr385 radical state (Intermediate II, E 5 ). This allowed us to analyse how the redox state of key catalytic components can influence the enzyme binding to AA.
The third COX 3 cycle, though not explicitly shown in the scheme in Figure 2, can be defined as a cycle formed by the enzyme states E 11 , E 19 , E 21 , and E 7 , involving reactions 47, 64, 33, and 53, respectively. The reactions in the COX 3 cycle are initiated in the result of AA binding to the enzyme being in a two-radical state [Fe(IV), PP*+, Tyr*] (state E 11 ). Note, that similar COX cycle was considered in the extended BC mechanism model [10,18], but its contribution to PGHS-1 functioning was not indicated.
Using similar approach, we introduced into our model the POX 3 and POX 4 cycles, intermediates of which differ from POX 1 intermediates by the occupancy of the COX binding site-it contains AA radical in the POX 3 cycle and PGG 2 radical in the POX 4 cycle. The POX 5 cycle was included to account for remaining POX activity, when the COX site is inactivated (see more detailed description below).
In our model, we also included an alternative route for initiation of COX reactions, which was proposed in an early version of the TC mechanism of PGHS-1 function [5]-a direct hydrogen atom abstraction from AA by protoporphyrin radical cation PP*+ and generation of the AA radical. In our scheme, this mechanism is described by the reaction of AA binding to the initial state E 1 resulting in the production of enzyme-substrate complex E 14 , [Fe(III), PP, Tyr, AA], which further converts into intermediate, [Fe(IV), PP, Tyr, AA*] (E 13 ), containing the arachidonyl radical AA* (reaction 56). Consideration of two alternative mechanisms of AA radical formation in our model allowed us to analyse their possible contributions to overall COX activity.

Modelling of PGHS-1 Self-Inactivation
PGHS-1 activity undergoes irreversible self-inactivation (SI) during catalysis [19]. The exact mechanism of this process is still not completely understood [3]. However, there is a consensus that SI most likely originates from the enzyme damage caused by oxidative reactions, coming from certain radical intermediates generated in PGHS-1 [20]. Potentially, any enzyme intermediate containing oxidising equivalents or radicals can cause the enzyme damage. One of the key questions discussed in the literature is which intermediate states can serve as the origin for POX and COX inactivation.
Based on the experimental data, we considered several possible mechanisms for enzyme SI in our model. These mechanisms included separate inactivation of the COX and POX activities, originating from enzyme intermediates with different redox state of heme group and containing one or two radicals in the catalytic domains [20]. All the SI processes accounted for in our model can be conditionally classified into the two following groups.
First, inactivation of POX activity resulting in fully inactivation enzyme (FIE state in Figure 2) was taken into consideration. According to the BC mechanism, SI of POX activity leads to the subsequent loss of COX activity, as POX reactions are required to regenerate Tyr385 radical state. In our model, this SI mechanism is described by reactions 43-46 and 54. Due to simplicity reasons, we did not Pharmaceuticals 2020, 13, 265 7 of 21 distinguish two stages in this process and approximated it by the following one-stage irreversible first order reactions: (a) Reactions 38-42, 53, and 67 originate from Intermediate I and lead to POX and further COX inactivation due to radical state of heme prosthetic group ([Fe(IV), PP*+]) [21,22]. In our model, this SI process can be initiated from seven catalytic states containing radicals such as E 2 , E 7 , E 11 , E 14 , E 19 , E 21 , and E 23 ; (b) Reactions 43-46 originate from Intermediate II (states E 5 , E 9 , E 15 , and E 20 ) and lead to POX and further COX inactivation due to the presence of oxidising equivalent Fe(IV) in heme.
Second, we considered the direct inactivation of COX activity with remaining POX activity which can be caused by the presence of radical states in the COX site in the model [19]. In our model, the remaining POX activity is described by the POX 5 cycle (reactions [25][26][27]. A similar POX cycle was considered in the model developed by Kulmacz et al. [18]. According to this hypothesis, we assumed that COX inactivation can originate from the enzyme intermediates containing either one or two radicals in the COX catalytic site. Reactions 58-62 and 68 originate from the states containing one Tyr385 radical in the COX site (E 5 , E 9 , E 11 , E 15 , E 19 , and E 20 ). Reactions 36 and 37 originate from the states with two radicals in the COX site, [Fe(IV),PP, Tyr*, AA*] (E 17 ) and [Fe(IV),PP, Tyr*, PGG 2 *] (E 4 ). It should be noted that the remaining POX activity can also further be inactivated due to damage caused by the presence of the radical state in heme (reaction 42).

Identification of a Set of Kinetic Parameters of PGHS-1
As a result of the joint fitting of three series of the experimental data (see Methods), a unified set of kinetic parameters was determined. The results of fitting are shown in Figure 3 and presented in Table 2. For the comparison of the obtained kinetic constants of PGHS-1 with the results of other authors, we also included experimental and theoretical estimations of the corresponding rate constants available from literature in Table 2.
Pharmaceuticals 2020, 13, x FOR PEER REVIEW 7 of 22 or two radicals in the catalytic domains [20]. All the SI processes accounted for in our model can be conditionally classified into the two following groups. First, inactivation of POX activity resulting in fully inactivation enzyme (FIE state in Figure 2) was taken into consideration. According to the BC mechanism, SI of POX activity leads to the subsequent loss of COX activity, as POX reactions are required to regenerate Tyr385 radical state. In our model, this SI mechanism is described by reactions 43-46 and 54. Due to simplicity reasons, we did not distinguish two stages in this process and approximated it by the following one-stage irreversible first order reactions: (a) Reactions 38-42, 53, and 67 originate from Intermediate I and lead to POX and further COX inactivation due to radical state of heme prosthetic group ([Fe(IV), PP*+]) [21,22]. In our model, this SI process can be initiated from seven catalytic states containing radicals such as E2, E7, E11, E14, E19, E21, and E23; (b) Reactions 43-46 originate from Intermediate II (states E5, E9, E15, and E20) and lead to POX and further COX inactivation due to the presence of oxidising equivalent Fe(IV) in heme.
Second, we considered the direct inactivation of COX activity with remaining POX activity which can be caused by the presence of radical states in the COX site in the model [19]. In our model, the remaining POX activity is described by the POX5 cycle (reactions [25][26][27]. A similar POX cycle was considered in the model developed by Kulmacz et al. [18]. According to this hypothesis, we assumed that COX inactivation can originate from the enzyme intermediates containing either one or two radicals in the COX catalytic site. Reactions 58-62 and 68 originate from the states containing one Tyr385 radical in the COX site (E5, E9, E11, E15, E19, and E20). Reactions 36 and 37 originate from the states with two radicals in the COX site, [Fe(IV),PP, Tyr*, AA*] (E17) and [Fe(IV),PP, Tyr*, PGG2*] (E4). It should be noted that the remaining POX activity can also further be inactivated due to damage caused by the presence of the radical state in heme (reaction 42).

Identification of a Set of Kinetic Parameters of PGHS-1
As a result of the joint fitting of three series of the experimental data (see Methods), a unified set of kinetic parameters was determined. The results of fitting are shown in Figure 3 and presented in Table 2. For the comparison of the obtained kinetic constants of PGHS-1 with the results of other authors, we also included experimental and theoretical estimations of the corresponding rate constants available from literature in Table 2. In the calibration procedure, we used two experimental datasets on PGG 2 and PGH 2 production by PGHS-1 treated by AA and phenol as RC [11] (Figure 3a-c). As a result of model fitting against the experimental data, the values of the dissociation constants of AA binding to the COX-site, K d,1 (reaction 1 in the COX 1 cycle) and K d,12 (reaction 12 in the COX 2 cycle) were estimated ( Table 2). The only experimental estimate available for the stage of AA binding to PGHS-1 is Michaelis constant K m ranged from 3 µM to 5 µM [1]. However, the Michaelis constant is an apparent value and generally depends on the concentration of the second PGHS-1 substrate RC (see discussion below) and therefore cannot be directly compared with K d,1 .
One of the key kinetic parameters of the PGHS-1 catalysis is the reaction rate of intramolecular electron transport [23] from Tyr385 to the heme group, leading to the formation of the Tyr385 radical. The obtained value for constant k 9 equal to 310 s −1 is in close agreement with theoretical estimates 300 s −1 [22] and 310 s −1 [4].
Fe(IV)PP* + Tyr → Fe(IV)PP + Tyr* k 9 310 s −1 300 s −1 [22]; 320 s −1 [4] Reactions in POX 2 cycle AA binding to the resting enzyme, E 1 POX reactions with adrenaline and H 2 O 2 Satisfactory quality of the fitting was achieved at the same values of kinetic parameters of COX activity for similar stages in all the three COX i cycles considered in the model. Note that the COX 1 , COX 2 , and COX 3 cycles differ only by the state of the heme prosthetic group in the POX site (see Section 2.1). The independence of kinetic constants of the COX reactions on the state of the heme group can indicate the absence of cooperativity between the POX and COX sites.
As a result of fitting, we estimated the rate constants for the different mechanisms of PGHS-1 self-inactivation ( Table 2). As discussed in Section 2.2, we considered three possible mechanisms of SI in the model. For SI originating from Intermediate I ([Fe(IV),PP*+]), we obtained the value of the rate constant k in1 = 0.6 s −1 , which is in good agreement with experimental data [4] (Table 2). For the reactions of COX inactivation (reactions 36, 37, 58-62, and 68) with the remaining POX activity (POX 5 cycle), we obtained the rate constant k in = 0.72 s −1 ( Table 2).
A small value of reaction rate k 8 of the reduction of protoporphyrin group radical state (Intermediate I, [Fe(IV), PP*+], E 2 ) by RC was obtained to be about zero k 8 (0.001 µM −1 s −1 ); that means that the Intermediate I is rather reduced by intramolecular process of electron abstraction from Tyr385 with radical formation than by exogenous RC.
To estimate uncertainty of the parameters obtained in the fitting procedure of the model against experimental data, we carried out global sensitivity analysis (GSA) of the model (see Methods). GSA revealed that the calibrated model showed uniform sensitivity to the most kinetic parameters ( Figure 4). The model is weakly sensitive to reaction rate k 12 of AA binding to the states in the POX 2 cycle (E 1 , E 2 and E 3 ), while it is sensitive to the dissociation constant K d,12 of the same reactions. The model also showed weak sensitivity to parameter k in2 , which corresponds to one of the SI reactions of the enzyme, but it is highly sensitive to parameters k in1 and k in characterising the other SI reactions considered in the model.

Contributions of the Different Enzyme Intermediates and Intra-Enzymatic Reactions to PGHS-1 Catalysis
The catalytic cycle of PGHS-1 ( Figure 2) represents a complex cycle network of intermolecular states with an intricate distribution of internal fluxes between them. To realise their contributions to the PGHS-1 catalysis, we carried out a comparative analysis of the population distribution of intermediate enzyme states and fluxes. This allowed us to estimate the relative impact of different elementary reactions and enzyme intermediates to the overall PGHS-1 function. We determined preferred ways of AA binding to the enzyme, identified key reactions producing intermediate PGG2 and final product PGH2, estimated the contribution of the different POX and COX cycles, and compared the role of different mechanisms of enzyme SI. The model also satisfactory described the third experimental dataset on a pure POX activity of PGHS-1 treated by H 2 O 2 and adrenalin as RC [12] (Figure 3d). As the result of fitting, another set of the kinetic parameters (k 5 , k 6 , and k 8 ) of the POX cycle was obtained, which differs from that identified for the POX reactions in PGHS-1 treated by AA and phenol as RC (Table 2). This result is consistent with experimental data on the pure POX activity of PGHS-1 in the presence of different peroxides and cosubstrates [5,6].
As seen from obtained agreement between theoretical and experimental data (Figure 3), the model with the unified set of kinetic parameters satisfactorily described the various experimental data such as kinetics of intermediate PGG 2 and final PGH 2 products, as well as consumption of substrate AA and cosubstrate RC. The model with the obtained set of parameters given in Table 2 was further validated and applied to analysis of regulatory mechanisms of PGHS-1 catalysis.

Contributions of the Different Enzyme Intermediates and Intra-Enzymatic Reactions to PGHS-1 Catalysis
The catalytic cycle of PGHS-1 ( Figure 2) represents a complex cycle network of intermolecular states with an intricate distribution of internal fluxes between them. To realise their contributions to the PGHS-1 catalysis, we carried out a comparative analysis of the population distribution of intermediate enzyme states and fluxes. This allowed us to estimate the relative impact of different elementary reactions and enzyme intermediates to the overall PGHS-1 function. We determined preferred ways of AA binding to the enzyme, identified key reactions producing intermediate PGG 2 and final product PGH 2 , estimated the contribution of the different POX and COX cycles, and compared the role of different mechanisms of enzyme SI.
As mentioned above, in the model, we considered two possible ways of initiation of the COX reaction-namely, two alternative ways of formation of the enzyme-substrate complex containing both Tyr385 radical and AA. The first corresponds to a generally accepted mechanism of the COX reaction initiation in the BC mechanism, when formation of Intermediate II state [Fe(IV), PP, Tyr*] (E 5 ) precedes the binding of AA (reaction 1). The second, alternative mechanism is the binding of AA to the COX site (reactions 30-32), which can be followed by the formation of the Tyr385 radical state E 9 through reaction 18 ( Figure 2).
The performed flux analysis showed that AA binds to the Intermediate II (E 5 ) state much more effectively than to the resting state (E 1 ), Intermediate I (E 2 ), and E 3 state. Maximum flux V 1 (t) exceeds maximum values of V 30 (t), V 31 (t), and V 32 (t) by an order of magnitude (see comparison of V 1 (t) and V 31 (t) in Figure 5a). The reaction rate of AA binding to the resting enzyme (V 31 ) has a high peak in an initial stage of the COX reaction, but its contribution to overall catalysis is negligible because of the short duration of the peak (1 ms). To determine what COXi cycles contribute significantly to the overall COX reaction, we compared the reaction rates of the PGG2 synthesis in various COXi cycles (Table 3) at the different concentrations of RC (10 μM, 100 μM, and 1000 μM). It turned out that the key COX cycle, in which the maximum amount of PGG2 is produced, is the COX2 cycle (reaction 8). In our model, the heme group of the COX2 cycle is in the ground ferric state [Fe(III), PP], as opposed to the kinetic models of PGHS-1 [5,9], where the heme group is assumed to remain in the unreduced ferril state [Fe(IV), PP] throughout the COX reaction (the COX1 cycle in our model).
Calculation of the reaction rates showed that the COX1 cycle plays a significant role only at low The results of modelling also showed that the rate constant of electron transfer from Tyr385 residue to heme group depends significantly on the presence of AA in the COX site (reactions 13 and 18). The rate constant k 9 of reaction 13 between enzyme states with the empty COX site equals to 310 s −1 , but the constant of reaction 18 between states with an occupied COX site k 11 is much lower (1.1 s −1 ).
Thus, reduction of the heme group due to Tyr385 radical formation in the presence of AA in the COX site is weaker than that in the case of the empty COX site. Therefore, substrate-enzyme complex [Fe(IV), PP, Tyr*, AA] (E 9 ) is formed mainly from Intermediate II according to the experimental data [5].
To determine what COX i cycles contribute significantly to the overall COX reaction, we compared the reaction rates of the PGG 2 synthesis in various COX i cycles (Table 3) at the different concentrations of RC (10 µM, 100 µM, and 1000 µM). It turned out that the key COX cycle, in which the maximum amount of PGG 2 is produced, is the COX 2 cycle (reaction 8). In our model, the heme group of the COX 2 cycle is in the ground ferric state [Fe(III), PP], as opposed to the kinetic models of PGHS-1 [5,9], where the heme group is assumed to remain in the unreduced ferril state [Fe(IV), PP] throughout the COX reaction (the COX 1 cycle in our model).
Calculation of the reaction rates showed that the COX 1 cycle plays a significant role only at low concentrations of RC < 10 µM (Table 3). This effect results from the increasing connection between COX 1 and COX 2 cycles that causes the reduction of heme group in the COX 1 cycle by RC (reactions 19, 22, 28, and 29) and a decrease of the role of the COX 1 cycle in the overall catalysis. Additionally, the COX 2 cycle significantly contributes to the activation of AA oxygenation reaction in comparison with that in the COX 1 cycle.
We also compared the rates of AA radical formation due to abstraction of a hydrogen atom from AA by the Tyr385 radical in both COX cycles (Figure 5b). The reaction rate V 6 in the COX 2 cycle was found to be significantly higher than the corresponding rate V 2 in the COX 1 cycle. We also compared the contributions of two different reactions of AA radical formation which can happen either due to the Tyr385 radical according to the BC mechanism of PGHS-1 catalysis or due to another oxidant-oxyferryl heme-according to the early version of the TC mechanism. The calculation showed that the reaction rate V 56 of AA* formation owing to oxyferryl heme electron-transfer was significantly lower than the corresponding reaction V 6 proceeding due to Tyr385 radical (BC mechanism) (Figure 5b). To analyse the main reactions contributing considerably to PGH 2 synthesis, we calculated the fluxes in the POX i cycles in which PGH 2 is produced ( Table 3). The main contribution to PGH 2 production turned out to be made by reaction 23 in the POX 4 cycle. In the PGHS-1 models developed on a basis of the BC mechanism [8], PGH 2 is produced in reaction 11 in the POX 1 cycle. By contrast, in our model, reaction rate V 11 is negligibly weak in comparison with reaction rate V 23 (Table 3). To analyse this effect, we calculated the population kinetics of the enzyme states. Figure 5c shows the time dependence of population of the states, maximum contributing to the PGHS-1 catalysis. These results showed that the POX 1 cycle contributes noticeably at the initial stage of PGH 2 synthesis, when the population of the resting state of the enzyme is maximum. In contrast, state E 12 in the POX 4 cycle works throughout all the catalysis, indicating that the POX 4 cycle is active up to full SI of PGHS-1. Moreover, state E 12 ([Fe(III), PP, Tyr, PGG*], being common for the COX 2 and POX 4 cycles, participates simultaneously in the production of the intermediate product PGG 2 (reaction 23) and the final product PGH 2 (reaction 8). As discussed above, these reactions make maximum contributions to the synthesis of PGG 2 and PGH 2 (Table 3). Thus, in our model, COX 2 and POX 4 cycles are the main cycles in PGHS-1 catalysis, while the COX 1 and POX 1 cycles contribute only at the initial stage of the PGHS-1 function.

Model Application to the Investigation of Regulatory Role of AA and RC in PGHS-1 Function
The calibrated model was applied to study the dependence of the PGHS-1 activity on substrate AA and cosubstrate RC concentrations. For comparison of the theoretical results with the experiments, we used the validation experimental data which were not used in the fitting procedure. The comparison of model and experiment results made without any optional fitting was a validation test for the predictive power of the model. This allows the model to be applied to the reliable analysis of in vivo experimental data on the NSAID action.
First, we applied the model to calculate the dependence of the oxygen consumption rate on AA concentration V O2 (AA). For this, we first calculated the kinetics of the oxygen consumption rate V O2 (t) and compared the results with the experimental data [8] (Figure 6a). As seen, the model satisfactory describes the experimental data on V O2 (t) and reproduces the time of the enzyme activity (about 20 s), which is in agreement with the SI time of PGHS-1 [1]. To plot the dependence of V O2 (AA), we took the maximum values of the rates V O2 (t) calculated at different AA concentrations. The theoretical dependence of V O2 (AA) along with experimental data [11] are given in Figure 6b.  An increase of Km with increasing RC concentration is accompanied by a change in the curve shape of VO2(AA) from hyperbolic to sigmoidal at high RC concentrations (Figure 6b). To investigate this positive cooperativity in PGHS-1 function, we plotted the Eadie-Scatchard dependence of ratio VO2/AA/PGHS on ratio VO2/PGHS (Figure 6d) and provided the corresponding experimental data [11]. As seen, the theoretical downward-concave curve VO2/AA/PGHS vs. VO2/PGHS significantly differs from the linear Eadie-Scatchard plot corresponding to Michaelis kinetics ( / = − ⁄ + ⁄ ) that confirmed positive cooperativity in COX activity of PGHS-1 [26].
Positive cooperativity in the COX reactions defines the activation threshold of the PGHS-1 catalysis which depends on both AA and RC concentrations [11,13]. Figure 7 shows the activation threshold curve on a plan of AA and RC concentrations which divides the range of AA and RC concentrations into two regions, i.e., PGHS-1 is inactive and active at values of AA and RC We also calculated V O2 (AA) at different RC concentrations (Figure 6b). As seen, RC inhibits COX activity of the enzyme, i.e., an increase in RC concentration leads to an increase of apparent Michaelis constant K m and decrease of the maximum oxygen consumption rate. This means that AA and RC compete each other for the COX site according to reactions 1 and 9 in the catalytic cycle ( Figure 2). Based on these results, we calculated the dependence of the Michaelis constant K m for AA on RC concentrations, K m (RC) (Figure 6c). The linear dependence K m (RC) looks like the dependence of the apparent Michaelis constant for the substrate on competitive inhibitor concentration, where RC plays a role of an inhibitor. The range of theoretical values of K m embraces a wide range of the experimental values of Michaelis constant, 3 µM [1], 5 µM [9], and 10.2 ± 1.5 µM [26]. Note that dissociation constants of AA binding with the enzyme obtained in our calculation is significantly lower (K d1,AA = 0.1 µM) than K m ( Table 2).
In the model, RC competes with AA for binding to states E 5 and E 15 in the COX 1 cycle (reactions 1 and 9) and COX 2 cycle (reactions 5 and 50), respectively. According with the BC mechanism, redaction of Tyr385 radical by RC leads to the inhibition of COX activity.
Moreover, the discussed mechanism may explain the distinction of the obtained kinetic constants of the reduction rate of Tyr385 radical depending on a state of the COX site. In the parameter set (Table 2), the rate constant k 10 of reactions 14 and 51 of Tyr385 reduction in states E 9 (COX 1 cycle) and E 20 (COX 2 cycle) with AA in the COX site equals approximately zero in contrast to the corresponding constant k 5 in reaction 9 and 50 of Tyr385* reduction in states E 5 (COX 1 cycle) and E 15 (COX 2 cycle) with the empty COX site ( Table 2). AA in the COX site is proposed not to permit RC to react with Tyr385*. For similar reaction 49 of Tyr385* reduction in states E 4 (POX 2 -cycle) with PGG2* in the COX site, the constant equals k 5 , which may mean that RC can reduce two radicals state in the COX site.
An increase of K m with increasing RC concentration is accompanied by a change in the curve shape of V O2 (AA) from hyperbolic to sigmoidal at high RC concentrations (Figure 6b). To investigate this positive cooperativity in PGHS-1 function, we plotted the Eadie-Scatchard dependence of ratio V O2 /AA/PGHS on ratio V O2 /PGHS (Figure 6d) and provided the corresponding experimental data [11]. As seen, the theoretical downward-concave curve V O2 /AA/PGHS vs. V O2 /PGHS significantly differs from the linear Eadie-Scatchard plot corresponding to Michaelis kinetics (V O2 /AA = −V O2 /K m + V max /K m ) that confirmed positive cooperativity in COX activity of PGHS-1 [26].
Positive cooperativity in the COX reactions defines the activation threshold of the PGHS-1 catalysis which depends on both AA and RC concentrations [11,13]. Figure 7 shows the activation threshold curve on a plan of AA and RC concentrations which divides the range of AA and RC concentrations into two regions, i.e., PGHS-1 is inactive and active at values of AA and RC concentrations lying below and upper the curve, respectively. The increasing activation threshold at increasing RC concentrations is determined by the inhibition of the autocatalytic COX reaction by RC. The PGHS-1 kinetics at AA and RC concentrations far from the activation threshold are close to the classical Michaelis-Menten kinetics that are typical for in vitro experiments on PGHS-1 kinetics and experiments on its inhibition by NSAIDs. The nonlinear effect in PGHS-1 kinetics becomes more pronounced near the activation threshold, in the range of low AA concentrations (0.05-1 µM), which are likely to relate to cellular conditions [1]. Below this threshold, the PGHS-1 is in a latent state and can be activated by either increasing AA or decreasing RC concentrations, as shown by the arrows in Figure 7. pronounced near the activation threshold, in the range of low AA concentrations (0.05-1 μM), which are likely to relate to cellular conditions [1]. Below this threshold, the PGHS-1 is in a latent state and can be activated by either increasing AA or decreasing RC concentrations, as shown by the arrows in Figure 7.

Effect of Reducing Cosubstrate on COX Activity of PGHS-1
We also applied our model to describe the experimental data on the effect of RC on oxygen consumption by PGHS-1 [10]. The results of the calculation along with the experimental data on the oxygen consumption rate, VO2, depending on phenol concentration are shown in Figure 8a. As seen, RC does not significantly affect the oxygen consumption rate at RC concentration less than 1000 μM

Effect of Reducing Cosubstrate on COX Activity of PGHS-1
We also applied our model to describe the experimental data on the effect of RC on oxygen consumption by PGHS-1 [10]. The results of the calculation along with the experimental data on the oxygen consumption rate, V O2 , depending on phenol concentration are shown in Figure 8a. As seen, RC does not significantly affect the oxygen consumption rate at RC concentration less than 1000 µM and considerably inhibits COX activity at RC concentration higher than 1000 µM.

Effect of Reducing Cosubstrate on COX Activity of PGHS-1
We also applied our model to describe the experimental data on the effect of RC on oxygen consumption by PGHS-1 [10]. The results of the calculation along with the experimental data on the oxygen consumption rate, VO2, depending on phenol concentration are shown in Figure 8a. As seen, RC does not significantly affect the oxygen consumption rate at RC concentration less than 1000 μM and considerably inhibits COX activity at RC concentration higher than 1000 μM.  For the comparison of these theoretical results with experimental data, we used experimental data corrected for the oxygen electrode signal [10]. Note that these data showed a slight increase of V O2 , at approximately 20% above the control values in the comparison with uncorrected data [10] and data obtained by another group [24], which have the pick of V O2 at about 2-3 times the control value. Our results showing no stimulation effect of RC on the COX reaction are closer to the corrected experimental and theoretical data [10]. According to our model, RC has an indirect stimulation effect on oxygen consumption. Calculation of O 2 consumption kinetics showed its increase at growing RC concentration (Figure 8b). In the model, the enhancement of O 2 consumption at RC concentration up to 1000 µM results from the extension of the functioning time of the enzyme due to its protecting from SI by RC. Increasing of the

Analysis of the Connection between COX and POX Activities
To analyse how our model describes stoichiometry between cosubstrate oxidation and arachidonic acid oxygenation during PGHS-1 functioning, we calculated the dependence of the ratio between the rate of cosubstrate (phenol) and arachidonate consumption (RC/AA) on the ratio between their initial concentrations (RC 0 /AA 0 ). Comparison of the theoretical results with experimental data [8] are shown in Figure 8d.
The dependence of RC/AA shows the extent of the interconnection between the POX and COX cycles in the course of PGHS-1 function. The obtained ratio RC/AA ranges from 0.5 to 1.3, which is consistent with the experimental range (0.6-1.3) [8], when RC/AA is in the region of 0.5-5. The theoretical data agree with the prediction of the BC mechanism of PGHS-1 catalysis (0-2) [8,19]. In our calculation, the dependence of RC/AA on RC 0 /AA 0 is saturation dependence rather than linear dependence [8].
As AA consumption defines PGG 2 production, and RC consumption relates to PGH 2 synthesis, the ratio RC/AA determines the extent of the PGG 2 conversion to PGH 2 . Because complete conversion of PGG 2 to PGH 2 requires two molecules of RC (RC/AA = 2), the obtained value of the ratio RC/AA = 0.5-1.3 shows that POX cycles are slower than the COX ones and complete conversion of PGG 2 is not reached. Thus, according to the model, some fraction of the intermediate product, PGG 2 , remains in solution at the end of the catalytic cycle and is accumulated there. The accumulation of PGG 2 in solution is proposed to be important for activation of PGHS-1 at cell stimulation.

Discussion
The microscopic network model of PGHS-1 catalysis was developed as the further extension of a series of the models developed previously on the basis of the branched chain mechanism of PGHS-1 kinetics [8][9][10]. Our model presents a network of the interconnected five POX and three COX cycles proceeding independently of each other and including all possible microscopic states generated in the POX and COX sites and possible transitions among these multiple states during catalysis. Analysis of the cycle structure of this network model of PGHS-1 was carried out by Yordanov and Stelling [27] on the basis of the developed graph-based approach to the investigation of combinatorial complexity of enzymatic kinetics and signalling network dynamics.
Model calibration was carried out on the basis of the experimental datasets obtained at the homogeneous experimental conditions of PGHS-1 treated by AA and phenol as RC [11] that allowed the consistent description of a majority of the kinetic data known for PGHS-1 catalysis under wide experimental conditions within a unified set of kinetic parameters. The calibrated model satisfactorily reproduced key features of the complex PGHS-1 dynamics, i.e., the interaction of COX and POX activities, processes of enzyme self-inactivation, its autocatalytic mechanism, and the activation threshold of the enzyme. Correct description of these characteristics of PGHS-1 in the framework of self-consistent kinetic model with the unified set of kinetic parameters is critical for in silico screening and prediction of the NSAID inhibition effects in different experimental assays and targeted cells.
The developed model was also applied to the experimental conditions, where PGHS-1 was treated by other substrates, adrenaline as RC, and H 2 O 2 as a substate of the POX site [12]. As a result, a new subset of the parameters responsible for binding of these substrates to the enzyme was identified ( Table 2). This showed that, for the further application of the model to experimental conditions different from those used in this work (e.g., other cosubstrates), it is necessary to identify a subset of the kinetic parameters related to the binding of other substrates to PGHS-1 and fix the rest of the kinetics parameters characterising intra-enzymatic kinetics obtained in the fitting procedure (Table 2). Moreover, identification of the kinetic parameters of hydrogen peroxide reaction with PGHS-1 in the model (Table 2) allows the model to be applied to the investigation of the PGHS-1 kinetics and NSAID effectiveness at different hydrogen peroxide levels in various cells and organs [15].
To validate and test predictive power of the developed model, we applied it to the description of the regulatory properties of the PGHS-1 and compared the results with additional experimental data not used in the calibration procedure. A study of the enzyme activity showed that AA and RC play various regulatory roles. We elucidated different mechanisms of inhibition and activation effects of RC on the enzyme. AA along with RC were shown to play an activation role, switching the enzyme from an inactive state to the functioning one. In this case, RC was found to influence the activation threshold of the enzyme. Note that detailed understanding and description of these regulatory properties are important for the correct analysis of platelets activation as well as NSAID effects on this process.
The developed model was applied to study intermediate enzyme states and reactions which contribute significantly to the synthesis of PGG 2 and PGH 2 . Taking into account two basic mechanisms of PGHS-1 function in the model, i.e., the branched chain [5,8] and tightly coupled [6,7] mechanisms, allowed us to estimate the contribution of each mechanism to the prostaglandin H2 biosynthesis.
The model supported the BC mechanism of PGHS-1 function, which assumes that Tyr385 radical in the COX site is a key component in interconnection between the COX and POX activities [8]. We also concluded that Tyr385 radical formation is required for effective binding of AA to the COX site. For AA to bind effectively, PGHS-1 first should be activated by a formation of Tyr385 radical due to POX activity of the enzyme. This effect is likely to result from a possible conformation change in the COX site by Tyr385 radical formation that can lead to an increase of the COX site affinity to AA. Conversely, AA in the COX site can suppress Tyr385 radical formation. Note that the experimental data indicates the possible mutual interaction between molecular compound in the COX site and Tyr385 radical formation. For example, experimental studies strongly indicate that the tyrosyl radical in inhibitor-treated PGHS-1 is probably not located at Tyr385 [28], and the alternative tyrosyl radical is not competent for COX catalysis. This data can serve an example of the effect of inhibitor molecules on the enzyme structure and internal dynamics particularly on the suppression of Tyr385 radical formation in the COX site [29].
As the result of the model fitting against the experimental data, dissociation constants K d,1 of AA binding to the states containing Tyr385 radical (E 5 and E 15 in the COX 1 and COX 2 cycles, respectively) was lower than the corresponding constants K d, 12 of AA binding to the states not containing Tyr385 radical (E 1 , E 2 , and E 3 states in the POX 1 cycle, see Table 2). A negligible contribution of reactions 30-32 to the catalysis did not permit us to carefully study the difference in dissociation constants of AA binding to various enzyme states. Additional experimental and theoretical investigations are needed for elucidation of the Tyr385 radical influence on enzyme affinity to AA.
In the model, we also investigated different mechanisms of PGHS-1 irreversible self-inactivation (SI) which was suggested to be an additional regulatory mechanism preventing the excess generation of signalling molecules by the enzyme in cells [3]. The SI mechanism in the model included separate inactivation of the COX and POX activities resulting from the enzyme oxidative damage caused by radical intermediates generated in the POX and COX sites of PGHS-1 [20]. The model with the obtained reaction rates of SI satisfactorily described stopping kinetics of PGH 2 production and full inactivation of the enzyme occurring within approximately 20 min (Figure 3a,b and Figure 6a).
The flux analysis of the network model showed the basic role of the COX 2 cycle which includes the ferric state of the enzyme [Fe(III), PP] (E 15 ). This conclusion contradicts the models where only the ferril state [Fe(IV), PP] (E 5 ) was considered in the COX reactions [5,8]. To remove this contradiction, further investigation of the connectivity between the COX 1 and COX 2 cycles is needed.
In the model, we considered an alternative mechanism for the COX reaction initiation proposed in the framework of the early version of the TC mechanism [6]. According to this mechanism, a direct hydrogen atom abstraction from AA by the ferryloxo-protoporphyrin radical cation of the heme group (reaction 56 from the E 14 state in POX 4 cycle) occurs, when AA is in the COX site. The flux analysis allowed us to compare this mechanism of AA radical formation with the alternative mechanism through Tyr385 radical, proposed in framework of the BC mechanism [5]. The obtained values of maximum reaction rates V 2 in the COX 1 cycle, V 6 in the COX 1 cycle, V 64 from the E 19 state, and V 56 showed that the main mechanism of AA radical formation is the BC mechanism. Thereby, the first stage of AA oxygenation including AA radical formation occurs mainly due to Tyr385* reduction in accordance with the BC mechanism.
To investigate the inhibition role of RC, we calculated the dependence of the apparent Michaelis constant K m (RC) for AA on RC concentration and showed an increase of K m with increasing RC concentration. The obtained change of K m can explain a wide range of its variation (3-5 µM) obtained in different experiments [1]. It is to be noted that there is no reliable data on RC nature, and its concentration in different cells and in vitro measurements of the dose dependence for NSAIDs are usually carried out at high RC concentration (e.g., 1000 µM of phenol). A change in the apparent Michaelis constant K m for AA under RC concentration variation affects the values of IC 50 for drugs measured at different RC and peroxide concentrations [16]. We also obtained dissociation constant K d,1 of AA (0.1 µM), which is an important value to estimation NSAIDs and AA competition for the COX site and correct analysis of drugs' IC 50 s in various experimental conditions [15].
The inhibition effect of RC on PGHS-1 catalysis observed in the model at high concentration of RC (>1000 µM, Figure 8a) can be explained by RC competition for Tyr385 radical in the COX site. RC molecule (e.g., phenol) is likely to be able to penetrate the COX site and reduce Tyr385 radical that abrogates oxygenation reaction. The same mechanism was proposed for the interaction of nitric oxide with PGHS-1, i.e., NO radical can terminate the COX activity by its reaction with the Tyr385 radical, thus forming 3-nitrotyrosine [30]. A similar mechanism was suggested to explain an inhibition effect of the phenolic antioxidant acetaminophen acting as reducing cosubstrate for the peroxidase protoporphyrin and preventing tyrosine radical formation in the COX site [31].
According to the experimental data, RC also causes a stimulation effect on PGHS-1 catalysis at low RC concentration [24,31]. Our model partly explained the stimulation mechanism of RC (phenol, acetaminophen, and others) by preventing the enzyme from suicide inactivation and extending the catalytic time of PGHS-1 (Figure 3a,b and Figure 8b). Thus, our model is recommended to be used in analysis of the PGHS-1 kinetics and NSAID effects in the range of high concentration of RC (>200 µM phenol). The further extension of the current model is required to study the stimulation role of RC in PGHS-1 catalysis at low concentration of RC.
Analysis of the oxygen consumption dependence on substrate and cosubstrate concentrations showed its sigmoidal shape in the range of low AA and high RC concentrations (Figure 6b) that points at positive cooperativity and nonlinear mechanism in PGHS-1 functioning. Positive cooperativity in the enzyme kinetics was also confirmed by the calculation of the Eadie-Scatchard plot that significantly differs from linear dependence and has a downward-concave shape (Figure 6d). The molecular mechanic of the positive cooperativity of the PGHS-1 function is not entirely understood. In the case of the dimer structure of PGHS-1, a possible mechanism of cooperativity can be determined by either allosteric interaction between two subunits of the holoenzyme or cooperativity within each monomer [24,26]. In our model, we proposed the independent function of the PGHS-1 monomers; thus, the manifestation of cooperativity in the model does not connected with its dimer structure. An increase in the extent of cooperativity with increasing RC and decreasing AA concentrations suggested that the POX activity and the reducing cosubstrate concentration are determinant factors in positive cooperativity in the COX reactions. Cooperativity of the COX activity in the model is linked to the autocatalytic effect of AA because it is observed profoundly close to the activation threshold of PGHS-1. The mechanism of cooperativity is related to positive feedback via the POX cycle, which is activated at high RC concentrations. Note that the positive cooperativity and noticeable activation threshold are absent in the inducible PGHS-2 isoform, whose kinetics is closer to the classical Michaelis-Menten kinetic [13,26]. This difference in the initiation kinetics of two isoforms of PGHSs in the range of cellular submiromolar concentration of AA suggests a different control mechanism in prostaglandin biosynthesis by constitutive PGHS-1 and inducible PGHS-2. In the case of constitutive PGHS-1, the activation threshold being adjustable to cellular conditions ensures a fast activation of latent enzymes and their inactivation by triggering of AA concentration to cellular demands near the activation threshold. The lower activation threshold and the absence of positive cooperativity in the initiation kinetic of PGHS-2 guarantees a fast and strong response of the enzyme to the cellular stimulus [13]. It is proposed that the distinctive regulatory properties of two isoforms of PGHSs reflect their different physiological roles in the same cells expressing both PGHS-1 and PGHS-2 and consume the same substrate and cosubstrate.

Parametrisation of the PGHS-1 Model
In accordance with the developed scheme for the PGHS-1 catalytic cycle (Figure 2), we developed a mathematical model for PGHS-1 catalysis. The model includes 28 ordinary differential equations (ODEs) formulated for 24 enzyme catalytic states (E i ) and six metabolite concentrations (AA, O 2 , RC, OC, PGG 2 , and PGH 2 ). The whole system of ODEs and conservation laws are presented in the Supplementary Materials. To run numerical simulation and calibration of the model, we used the software package for kinetic modelling DBsolve 6.1 [32] and the Systems Biology Toolbox for MATLAB (The MathWorks, Inc.).
The reaction rate expressions for all the reactions were formulated according to the mass action law with the rate constants being the parameters of the model (see Supplementary Materials). To identify the model parameters, we fitted the solutions of the ODEs against the experimental data on enzyme kinetics of PGHS-1 available from the literature.
The model includes 66 reactions and 68 corresponding reaction rate constants. The experimental data published in the literature are not sufficient for the reliable identification of all the model parameters. The lack of the data is mainly due to some peculiarities of the experimental design generally used for kinetic studies of PGHS-1 catalysis. The majority of such experimental studies included either the time courses for metabolite concentrations or reaction rates of their consumption or production such as AA or O 2 consumption, RC oxidation, and PGH 2 production. Whereas these characteristics are relatively easy to study experimentally, the measurement of kinetic dependencies for individual enzyme intermediates (catalytic states E i in our model) requires the application of different experimental techniques, mainly EPR spectrometry for the radical intermediates detection [14]. To reliably identify 68 reaction rate constants, we would ideally need time series data for the majority of enzyme catalytic states considered in the model. Taking into account the lack of such data, we proposed the following way to reduce the number of parameters in the model.
In general, all the reactions in the developed catalytic cycle shown in Figure 2 can be classified into several groups, as some of them are of similar types. For example, we grouped together the reactions of AA binding to the COX site (reactions 1, 30, 31, and others) or reduction of PGG 2 to PGH 2 (reactions 11, 16, 20, and others). The reactions within the groups differ from each other by the enzyme catalytic state participating in corresponding reactions. For example, AA can bind to the enzyme in the resting state [Fe(III), PP, Tyr] (E 1 , reaction 31) and ferryl states [Fe(IV), PP, Tyr] and [Fe(IV), PP*+, Tyr] (E 3 and E 2 states in reactions 30 and 32, respectively). Alternatively, AA can bind to the enzyme state containing Tyr385 radical (E 5 , E 15 , and E 11 states in reactions 1, 5, and 47, respectively). Therefore, it was assumed that the reactions of a similar type would have similar rate constants.
However, in general case, kinetic constants of the elementary reactions taking place in one of the catalytic sites can depend on the state of all other catalytic components, located in the nearest proximity and involved in the further catalytic steps. For example, there is experimental evidence that AA binding to the COX site may depend on the presence/absence of the radical on Tyr385 [18]. PGHS-1 is a bifunctional enzyme, and although the reactions traditionally referred to as COX activity and POX activity are catalysed in different sites, there are indications that the processes in one site can influence the rate of the processes in another. This phenomenon is traditionally discussed in the context of the cooperativity of PGHS-1 [10].
We tried different approaches in grouping the model parameters and analysed the results of fitting of the corresponding models to the same set of experimental data. Analysis of the fitting results allowed us to hypothesise which reaction rates do not depend much on the processes taking place in the neighbour site and therefore which kinetic constants can be assumed to be equal to each other. This analysis also allowed us to make conclusions on the cooperativity between COX and POX activities. The approach we used allowed us to significantly reduce the number of parameters subject to fitting. The best fitting of the model to experimental data was achieved with the set of 18 parameters given in Table 2.
When identifying the model parameters, we were trying to find a unified set of parameters which would allow us to consistently describe a majority of the in vitro kinetic data known for PGHS-1 catalysis. Identification of the model parameters and their further validation was performed as follows. We chose several experimental data sets obtained in different research laboratories to fit the model parameters against the data and then validated the calibrated model by checking its applicability to the description of experimental data which were not taken in the fitting procedure.
To avoid ambiguity, for the initial stage of the fitting, among all kinetic information available for PGHS-1, we chose three homogeneous datasets which were obtained for the same object-purified PGHS-1 from sheep seminal vesicles-under approximately the same experimental conditions (temperature and pH). The first set included kinetic data on PGG 2 and PGH 2 production (points in Figure 3a,b) measured at the different concentration of RC, phenol [8]. As the second dataset, we used the kinetic data on AA consumption (points in Figure 3c), obtained for the different AA concentrations at the fixed concentration of the RC phenol [11]. The experimental conditions of the third dataset [12] differed from the first two. In this experiment, the kinetics of the separate POX activities of PGHS-1 treated by H 2 O 2 as the substrate of the POX site and adrenaline as RC were measured in the absence of substrate AA. These data were used to identify kinetic parameters k 5 , k 6 , and k 8 of the POX reactions for PGHS-1 treated by H 2 O 2 . The obtained values of these parameters were different from those identified from the experimental dataset on the kinetics of PGHS-1 treated by phenol as RC and AA (see Table 2). Note that the calibration of the model based on two different experimental conditions allows the model to be applied to the analysis of experimental data obtained for PGHS-1 treated by the different RCs and different substrates of the POX site (e.g., PGG 2 and H 2 O 2 ).

Global Sensitivity Analysis of the Model
Global sensitivity analysis (GSA) of the model to the kinetic parameters obtained in the fitting procedure was based on the Morris GSA method for the calculation of elementary effects [25,33]. The resulting Morris ranking is calculated as where µand σare two sensitivity measures in the Morris GSA method, i.e., µassesses the impact of a parameter on the model output, and σdescribes the non-linear effects and interactions of the parameters. As the model output in the calculation, we took an area under the curve (AUC) of PGH 2 time course.

Conclusions
The key aim of this work was the development of a reliable and validated kinetic model of PGHS-1 which reproduces the main features of its complex kinetics and can be used as a computational tool for accurate prediction of inhibition effect of NSAIDs and their combinations in a range of cellular conditions. The cycle-based network model integrates both the current models of the molecular mechanism of PGHS-1 catalysis and a large dataset of experimental data on the kinetics and regulatory properties of PGHS-1. This integrative approach helped us to overcome the limitation of the previously developed PGHS-1 models mainly applicable to the description of individual in vitro data in a specific range of substrate and cosubstrate concentrations. The developed model of PGHS-1 catalysis is a platform for the further development of a kinetic model of inducible isoform PGHS-2 (COX-2) which is expressed together with PGHS-1 in vascular endothelia cells and other organs and is a target of COX-2 selective NSAIDs. Development of the kinetics models of two isoforms within a unified integrative approach would allow for the investigation of the different physiological roles of these isoforms in the same cells and analysis of the effect of NSAID inhibition on the balance of pro-thrombotic (thromboxane) and anti-thrombotic (prostacyclin) prostaglandins in platelets and endothelial cells.