Stochastic Simulations of the Cytochrome P450 Catalytic Cycle - The

An internally developed C-language program of the next reaction method is used .... to the degree that the normal catalytic feature of P450 enzymes is...
0 downloads 0 Views 249KB Size
Download PDF
J. Phys. Chem. B 2007, 111, 4251-4260

4251

Stochastic Simulations of the Cytochrome P450 Catalytic Cycle Yonghua Wang,*,† Yan Li,‡ and Bin Wang† Key Laboratory of Mariculture and Biotechnology of the Ministry of Agriculture, Dalian Fisheries UniVersity, Dalian 116023, China, and School of Chemical Engineering, Dalian UniVersity of Technology, Dalian 116012, China ReceiVed: February 13, 2007; In Final Form: February 19, 2007

Cytochrome P450 enzymes involve a complex reaction cycle which has been described, for the first time, by a stochastic simulation of the system in the present work. A series of models are developed for a basic catalytic cycle, employing a set of microscopic rate constants for the oxidation of p-alkoxyacylanilides catalyzed by the cytochrome P450 1A2. By analyzing the effects of low concentrations of enzyme and substrate on the system, and the dependence of the system on several rate constants, it is discovered that the system evolves along relatively stable patterns from its initial state, as indicated from different runs of simulations. Strong fluctuations appear at the entrance and exit of the pathway, with very weak fluctuations in the middle sections of the cycle. Although noises are apparent when the reactant populations are very low, basically, the fundamental feature of the P450 cycle based on a microscopic view is that it is deterministic in nature. Meanwhile, the mathematical models we developed are qualitatively validated by a comparison with those experimental results of the P450 cycle. The findings of this work will be helpful for a further deeper understanding of the catalytic mechanism of cytochrome P450 enzymes.

1. Introduction Cytochrome P450 enzymes are major catalysts involved in the oxidation of the xenobiotic chemicals,1,2 a significant focus of scientists in the areas of toxicology, drug metabolism, pharmacology, and even theoretical chemistry. In the past decade, an impressive number of investigations, both experimental and computational, have been devoted to the better understanding of the P450-mediated catalytic mechanism.2,3 P450 catalysis involves a complex reaction pathway, probably including nine steps independent of protein conformational changes (Scheme 1).1 Elucidation of the reaction cycle is complicated by the difference of specific P450 enzymes in the influence of substrates, the sequence of events and the ratelimiting steps, as well as the organization of different isoforms in the microsomal membrane.2 The history of investigating the catalytic cycle goes back to shortly after the discovery of the P450 system.4-6 These efforts by various methods have continued to this day, with research of different P450s and reactions following multiple courses.7,8 An important attempt of the efforts might be the utilization of nonlinear ordinary differential equations (ODE) to establish a mathematical model based on the reaction kinetics for P450s.9,10 This approach can provide us with mathematically well-founded and tractable interpretations regarding pathways by the enzyme reactions. More recently, we also proposed a system-theoretical approach, a time-dependent metabolic control analysis, to the analysis and quantitative modeling of the CYP450 catalytic pathway.11 This provides a theoretical enlightenment for us to assess the transient response of the system to perturbations. This biochemical modeling and simulation largely relies on the deterministic methods describing the evolution of the popula* To whom correspondence should be addressed. E-mail: yhwang@ dlfu.edu.cn. † Dalian Fisheries University. ‡ Dalian University of Technology.

SCHEME 1: The P450 Catalytic Cycle for the Stochastic Modelinga

a E, P450; S, substrate; ESi (i ) 1, 2, ... , 5), intermediate complex; EP, Fe3+-ROH complex; P, product.

tions with differential equations. However, this law is valid only if the behavior of a particular reactant is ignored and only the average behavior of the populations is considered. This ideal situation is certainly rarely matched in cell biology, where usually a small number of reacting entities (in dilute solution) causes the individual behaviors to be significant. Recently, many studies have reported the occurrence of stochastic fluctuations and noise in living systems.12,13 Intrinsic stochasticity is inherent to the system, arising due to very low

10.1021/jp071222n CCC: $37.00 © 2007 American Chemical Society Published on Web 04/03/2007

4252 J. Phys. Chem. B, Vol. 111, No. 16, 2007 reactant populations. However, extrinsic stochasticity also originates, due to the random variation of one or more environmental factors such as the temperature and concentrations of the reactant species. Stochastic fluctuations have been found to play important roles in many biologically important reactions.14,15 Following this finding, considerations of the occurrence of stochasticity or randomness during the modeling of the cellular pathways have also been focused on recently.12 A system with stochastic characteristics is in favor of biological simulation using stochastic algorithms. However, deterministic approaches, like the ODE-based methods, take a macroscopic view of the biological reactions, describing averaged dynamics, since they assume that the molecules are uniformly distributed. Thus, it is not possible to capture emergent phenomena that arise from inherent randomness. In contrast to those concentrations that are real numbers in the deterministic approach, stochastic simulations which employ concentrations of discrete numbers take a mesoscopic view of a system and keep track of every molecule in the system, which is close to the real world accuracy. Although a large amount of data focusing on the P450 topics have been reported, these reports mainly rely on deterministic approaches and do not provide information on the P450 mechanism underlying stochasticity. In this work, a stochastic version of the P450 catalytic cycle is developed, so as to investigate if certain stochasticity exists in the system, and, if so, how it affects the behavior of the enzyme. The study is broadly divided into two sections, with the first section introducing the Gillespie algorithm and the second one presenting the stochastic results of the P450 cycle, and discussing the effect of stochasticity on the kinetics of the cycle. 2. Methods Gillespie’s Algorithm. The intrinsic and extrinsic stochasticity of the biological system lead to the development of stochastic modeling techniques. Two widely used techniques for the stochastic simulations are Monte Carlo simulations16 and Gillespie’s method,17 with the latter being the most important stochastic algorithm. The advantage of Gillespie’s algorithm is that it generates an ensemble of trajectories with correct statistics for a set of biochemical reactions. In contrast to the deterministic formulation of chemical kinetics, the algorithm remains exact for arbitrary low numbers of molecules. In 2000, Gibson and Bruck proposed a method, the next reaction method (NRM), as an enhancement of Gillespie’s first method.18 The NRM is efficient to improve the performance of the Gillespie algorithms substantially while maintaining exactness of the algorithm. Gillespie algorithms require enormous computational time for a system with a large number of reactions. The next reaction method uses only one random number per iteration as compared to the M random numbers in the first reaction method. The reduction in generation of random number and better bookkeeping of the data has improved the performance of Gillespie algorithms significantly. Gillespie algorithms have been successfully applied to in silico biological simulations recently,19,20 which will also be used for the analysis of the P450 catalytic mechanism in this work. An internally developed C-language program of the next reaction method is used in the present study, using the following implementations: 1. Initialize • Set the initial number of species, t ) 0, and generate a dependency graph G which was built on the basis of Table 1. • Calculate the propensity, fj(x), for all j. • For each j, generate a putative time, τj, according to the exponential distribution with parameter fj(x).

Wang et al. • Store τj in queue P. 2. Let u be the reaction whose putative time, τu, stored in P is least. Set τ ) τu. 3. Update the states of the species to reflect the execution of reaction u. Set t ) τ. 4. For each (R, υ) in G, • update fR • if R * υ, set τR ) fR,old/fR,new(τR - τ) + t • if R ) υ, generate a random number, r, and compute τR according to the following equation: τR ) (1/fR(x)) log(1/r) + t • replace the old τR in P with the new value Go to step 2 for a new iteration. The Stochastic Model of the P450 Catalytic Cycle. An abbreviated scheme of the P450 pathway according to Guengerich9 is modeled (Scheme 1), with a total of 9 types of biochemical reactions as well as 11 species involved in the cycle. In this work, the initial numbers of P450 enzyme and its substrate, including the kinetic constants, are assigned according to refs 9 and 21 or theoretical deviations. The rates of the individual steps of catalytic cycle of rat P450 1A2 for phenacetin oxidation are the following: k1 ) 6000, k-1 ) 120 000, k2 ) 700, k3 ) 6000, k-3 ) 6, k4 ) 700, k5 ) 115, k6 ) 30, k7 ) 660, k-7 ) 6000, k8 ) 50, k9 ) 50. All rates are in min-1 (first order) or µM-1 min-1 (second order). In the stochastic modeling, all of these kinetic parameters have been expressed in terms of elementary first- and second-order reactions. Table 1 shows the stochastic version of the P450 cycle. The first column is the reaction number. The second one is the reaction sequences. The propensity of each reaction is given in the third column. The propensity is a product of the reaction rate and a function of the number of the molecules involved. 3. Results and Discussion In this work, we applied the stochastic algorithm for simulation of P450 catalytic mechanisms. Our attempt for this stochastic description of the P450 cycle is in line with intuitions: (1) The numbers of individual P450 enzyme and substrate involved in reaction per cell are small. Stochasticity becomes pronounced when the number of reacting species is low. It is known that some P450 isoforms such as the cytochrome P450s 1A2, 3A5, and 2C1822 occur in low numbers in hepatic cells, and generally the xenobiotic chemicals (substrates) are also low in ViVo. Particularly, the P450 large catalytic pocket could only accommodate one or two substrates.23 (2) Limited diffusion effects are apparent due to the structural organization of the cytoplasm and the macromolecular crowding of P450s (heterogeneous).24 The heterogeneous intracellular environment makes the stochasticity inevitable. It has been found that the time evolution of such a chemical system with stochasticity is a Markov type of stochastic process.25 However, in a deterministic simulation space, it is not possible to capture emergent phenomena that arise from inherent randomness. All of these allow us to perform a stochastic simulation of the P450 catalytic cycle. Temporal Trajectory. The initial concentrations of enzyme and its substrate are [E]0 ) 30 nM (estimated)22 and [S]0 ) 30 nM, respectively. In the modeling, a series of realizations were carried out and averaged before the property descriptions of the stochastic system were obtained, with a purpose of illustrating the overall trend of the species, which is free of the variability in single realizations. Figure 1A illustrates a temporal trajectory of all species of the cycle, and Figure 1B shows their distributions in all time ranges. One can find that the P450 enzyme (E) changes insignificantly in a long time run and thus it could be assumed constant, and the substrate (S) attenuates

Cytochrome P450 Catalytic Cycle

J. Phys. Chem. B, Vol. 111, No. 16, 2007 4253

TABLE 1: The Stochastic Model of the Cytochrome P450 Cyclea reaction no.

reaction k1

propensity of reaction

depends on

P1 ) k1 × (#E) × (#S)

E, S

k-1

P2 ) k-1 × (#ES)

ES

k2

P3 ) k2 × (#ES1)

ES1

P4 ) k3 × (#ES2) × (#O2)

ES2, O2

k-3

P5 ) k-3 × (#ES3)

ES3

k4

P6 ) k4 × (#ES3)

ES3

k5

P7 ) k5 × (#ES4)

ES4

k6

P8 ) k6 × (#ES5)

ES5

P9 ) k7 × (#EP)

EP

P10 ) k-7 × (#P) × (#E)

P, E

k8

P11 ) k8 × (#ES4)

ES4

k9

P12 ) k9 × (#ES5)

ES5

P13 ) k10 × (#ES4) × (#H+)

ES4, H+

k11

P14 ) k11 × (#EC) × (#H+)

EC, H+

k12

P15 ) k12 × (#EC) × (#H+)

EC, H+

1

E + S 98 ES

1′

ES 98 E + S

2

ES1 98 ES2

3

ES2 + O2 98 ES3

3′

ES3 98 ES2 + O2

4

ES3 98 ES4

5

ES4 98 ES5

6

ES5 98 EP

7

EP 98 P + E

7′

P + E 98 EP

8

ES4 98 ES1 + H2O2

9

ES5 98 ES1 + H2O

10*

ES4 + H + 98 EC

11*

EC + H + 98 ES5 + H2O

12*

EC + H + 98 ES1 + H2O2

k3

k7

k-7

k10

a # Q denotes the number of Q. Steps 10, 11, and 12 are three reactions newly added in Scheme 1. EC: Fe3+-OOH-‚RH. Step 3, O binding, 2 is generally considered to have a low Kd, as judged by the ready oxidation of ferrous P450.34 This step is a step assumed to be determined by diffusion, which is treated as being near diffusion-limited in kinetic analysis. We assume the concentration of O2 in every reaction of simulation was 1 nM.

gradually to extinct. These results conform to the experimental observations, since the enzyme concentration will recover to its initial state after a reaction cycle is performed. By an observation of the overall trend of the system, an interesting phenomenon was discovered that strong fluctuations appear at the entrance and exit of the pathway, with very weak fluctuations occurring in the middle sections of the cycle. The fluctuations for both E and EP (complex) are clearly evident; instead, the product (P) increases steadily and its fluctuation can seldom be observed. In addition, ES1 fluctuates slightly, whereas ES2 bursts with strong oscillations at the beginning of its formation and then attenuates. The middle section reveals very stable patterns for ES3, ES4, and ES5. These observations might indicate that the degree of the fluctuations is strongly dependent on the positions where they are present in the cycle, which may be a new catalytic mechanism of the P450 cycle. Intuitively, strong fluctuations appearing at the entrance of the cycle are vital, as they may play a causal role in the initiation of the substrate seizure for P450 enzyme. However, the weak fluctuations in the middle section are also important, which may be closely related to the stability of the P450s’ functioning. When modeling with different numbers of initial S and E, and for each specific number of molecules, it was discovered that the trend of the cycle obtained from the average results over several runs is almost identical with the ones of the single run. In principle, when the number of molecules becomes small,

the variability of the molecular populations in the biological reactions increases. To identify whether there is observed variability, the initial concentrations were assigned as [S]0 ) 5 nM and [E]0 )10 nM (nearly 10 molecules/cell), respectively. This simulation yields the same average behavior with no qualitative difference compared with the results of Figure 1. Steady changing patterns are found for S, H2O, and H2O2, yet for E, obvious oscillations are observed. As a matter of fact, the current concentrations of reaction species have been assigned very low. However, it is discovered that the fluctuations are less prominent, and for most species, they even tend toward the deterministic solution (Figure 2). As compared with the situations with lower concentrations (Figure 2), fluctuations become weaker when the reaction species are in higher concentrations (Figure 3, [S]0 ) 150 nM and [E]0 ) 150 nM). Considering that in hepatic cells the concentrations of P450s like P450 3A426 and its substrate are both high, the metabolic reactions containing fluctuations with uncertainty could be weak. From a comparison of Figures 1, 2, and 3, no significant difference in the temporal evolution patterns of all species is observed. Although fluctuations to some extent exist in the system, in a long time run, all of the temporal evolution patterns tend to be deterministic, and the system is stable for different initial conditions. Meanwhile, no periodic oscillations like those occuring in many other cellular pathways (cell cycle, etc.) are observed in the simulations of this cycle (Figures 1-3).

4254 J. Phys. Chem. B, Vol. 111, No. 16, 2007

Wang et al.

Figure 1. (A) Stochastic simulation using Gillespie’s next reaction method for the P450 catalytic cycle. The kinetic constants of the rat P450 1A2 for phenacetin oxidation are the following: k1 ) 6000, k-1 ) 120 000, k2 ) 700, k3 ) 6000, k-3 ) 6, k4 ) 700, k5 ) 115, k6 ) 30, k7 ) 660, k-7 ) 6000, k8 ) 50, and k9 ) 50. All rates are in min-1 (first order) or µM-1 min-1 (second order). Both the initial enzyme (E) and substrate (S) concentrations were 30 nM. (B) The distribution graph of each molecule in the P450 cycle in the whole time range. X axes represent the number (concentration) of a certain molecule, and Y axes represent the ratio of the distribution of each concentration in the total concentration of a certain molecule.

Therefore, the stochastic simulations reveal that the P450 cycle system is inherently deterministic, and not random. Cytochrome P450 enzymes are found throughout nature, from bacteria to human. The function of the cytochrome P450

superfamily is to facilitate the elimination of xenobiotic and endogenous compounds by metabolizing them in ViVo.27 Therefore, a deterministic reaction way should be a prerequisite for the normal protection function of the enzyme for various

Cytochrome P450 Catalytic Cycle

J. Phys. Chem. B, Vol. 111, No. 16, 2007 4255

Figure 2. Stochastic simulations using Gillespie’s next reaction method for the P450 catalytic cycle with the initial concentrations of [S]0 and [E]0 assigned as 5 and 10 nM, respectively, and all other conditions assigned the same as those in Figure 1.

Figure 3. Stochastic simulations using Gillespie’s next reaction method for the P450 catalytic cycle, with both the initial enzyme (E) and substrate (S) concentrations being 150 nM. Other conditions are the same as those in Figure 1.

organisms. It is revealed that the metabolic function of the enzyme is in a deterministic way and no periodic fluctuations between high and low numbers to the system occur. Thus, in this way, its substrates can be removed by being oxidized eventually (Figures 1-3). Under low concentration conditions, simulations reveal that noise or fluctuations may have impacts,

but not to the degree that the normal catalytic feature of P450 enzymes is changed. Simulations of the Formation of H2O2 and H2O. Step 8 represents the conversion of (FeO2)+ to H2O2 (alternatively, H2O2 could be formed from Fe2+O2; such a possibility has not been included). The formation of H2O is presumed to occur

4256 J. Phys. Chem. B, Vol. 111, No. 16, 2007

Wang et al.

Figure 4. Effects of different k6 values (min-1) on the production patterns of the final product of the P450 cycle in time evolution.

only from a reduction of the species labeled (FeO)3+ in Scheme 1. Accounting for that, the formation of H2O2 and H2O is critical in this system; although the present simulations involve relatively simplified mechanistic considerations to the formations of the two molecules, they are still of value. In addition, the deficiency of most current work with the P450s is their lack of ability to measure the H2O2 and H2O directly; theoretical attempts are possibly the major way to understand the behaviors of the two molecules. In this modeling, steady changing patterns with very weak variability of the formation of H2O2 and H2O are observed, as shown in Figure 1. Their physical interpretations would be that the access of water molecules into the pocket needs to be well controlled. On one hand, water molecules are important for protonation, but on the other hand, too much water leads to uncoupling. Strong fluctuations of H2O molecules would be unreasonable for an efficiency of the catalytic system. By the way, the modeling of H2O2 and H2O also presents a comparison of the stochastic approach with the deterministic one for the P450 system. For stochastic modeling, it is feasible to take a mesoscopic view of the system which can keep track of every molecule in the system. As to a discrete number of molecules, the concentration change corresponds to certain reaction events. The algorithm determines the exact time when the individual molecular reactive encounters occur. Here, we present a snapshot image of the changing patterns of species in the cycle in time evolution. For H2O2, there are two graphs concerned with its formation (Figure 1), where the graph with time varying from 0 to 15 000 ms is a snapshot of the one with time varying from 0 to about 4 × 105 ms. Following the time pathway of the molecule, one can discover that in every time step only one H2O2 molecule (1 nM ≈ 1/cell) is generated, and every other species of the cycle is exactly in the same situation. In this way, the stochastic modeling presented detailed information of the time courses of all species. However, the deterministic approach cannot offer such information, as its basic assumption is that the concentrations of species vary continuously over time. In other words, the concentrations are of the real numbers instead of the discrete numbers in the stochastic

simulations, and also the variation of concentrations shows no discontinuities or “jumps”. However, even with nine steps, the paradigm presented in Scheme 1 is an overall simplification, with some important (stable) intermediates possibly being missed,28 as the rate constant for each step of the cycle is still a hard problem to solve under the current experimental conditions. Therefore, further stochastic studies of the cycle based on a more extensive experimental study of certain key rate constants are still needed. Kinetic Parameter Effects. To investigate the dynamic effects of kinetic parameters on the cycle, some additional computational analyses were performed. Previous research revealed that a change in k6 (which was considered as the most sensitive to changes in terms of losing the curve fits) produced dramatic effects on the cycle.10 Therefore, in our work, several changes have been made for k6, in turn, that is, k6 ) 6, 60, or 1000 min-1, respectively, to investigate its effects on the formation of product. Figure 4 illustrates the analysis of k6. When the reaction rate increases, the product concentration also rapidly increases (Figure 4). It seems sensible that the larger the rate constant is, the more probabilities the collision of the reaction species for the product formation have. These computational analysis results are consistent with those observed in the experimental studies. In addition, analyses were also conducted for other related parameters, that is, k4 and k5. They were changed synchronously (k4 ) 350 min-1 and k5 ) 90 min-1, smaller than their previous data in Figure 1), and other conditions were kept the same as those in Figure 1. The reason why we only changed these two parameters is due to the fact that these changes are experimentally feasible and can satisfy several criteria, including the isotope effects, the formation rates of the reduced O2 species, as well as the V-vs-S plots.10 The model describes a system in which the changes in k4 and k5 can affect the dynamic behaviors of the system. However, as seen from Figure 5, there is no significant difference between this figure and Figure 1, with only a relatively longer clearance of substrate observed in the former. We know that for the P450 cycle the magnitude of the reaction rates of its individual steps

Cytochrome P450 Catalytic Cycle

J. Phys. Chem. B, Vol. 111, No. 16, 2007 4257

Figure 5. Stochastic simulations using Gillespie’s next reaction method for the P450 catalytic cycle, with all conditions the same as those in Figure 1, except that k4 ) 350 min-1 and k5 ) 90 min-1.

Figure 6. Stochastic simulations with k-1 assigned as 12 000 µM-1 min-1 for parts A and B and 1200 µM-1 min-1 for parts C and D.

affect the turnover efficiency of the enzyme, but they may not affect the whole evolution trend of the cycle, as shown in Figures 1 and 5. The ratio value of k-1/k1 (the substrate binding) used in these models is similar to the Kd value estimated by the fluorescence titrations with substrates (20 µM).9 The k-1 reaction seems to be the fastest in the model and intuitively implies that it is difficult to get the catalytic cycle starting. Therefore, we tried

different values for k-1, that is, 12 000 and 1200 min-1 µM-1, to see its effects on the whole cycle. For simplicity, Figure 6 only shows the results of E and S. The results are similar to those shown in Figure 1, except that the larger k-1 is, the more slowly the enzyme clears the substrates (Figures 1 and 6). This is explainable, since the first reaction is a going in and off one, which may be affected by the diffusion-limited encounters of substrates with ligands. We have not made an attempt to change

4258 J. Phys. Chem. B, Vol. 111, No. 16, 2007

Wang et al.

Figure 7. Stochastic simulations using Gillespie’s next reaction method for the P450 catalytic cycle, with the introduction of the intermediate Fe3+-OOH- (EC) to Scheme 1. All three second-order reactions were set at 6000 µM-1 min-1, with pH 7 ([H+] ) 100 nM). Other conditions are the same as those in Figure 1.

the value of k-1 (the off rate), which always takes on different rates seen in different P450 systems,29,30 and this, at the same time, suggests the complexity of the systems studied. Although a few general catalytic mechanisms appear to be operative for most of the reactions catalyzed by the P450s,31 certain features such as the rate-limiting steps and the substrate interactions may vary considerably depending on the particular P450 and reactions involved.32 From the above simulation, we find that, to some extent, the dynamic time-evolution patterns for the specific reaction p-alkoxyacylanilides catalyzed by cytochrome P450 1A2 are affected by the reaction steps. However, it is still unclear how other P450 isoforms behave with respect to the changes of parameter values adopted in building the stochastic model. Our present analysis of the effects of the individual rate constants provides a framework for consideration of other P450 reactions. Stochastic Bifurcation Prediction. The P450 catalytic pathway clearly exhibits a nonlinear behavior because of the complex underlying mechanism of the interactions, such as the enzymatic actions and feedback loops (steps 8 and 9, the H2O and H2O2 formation). Due to the nonlinearity as well as the low concentration conditions, the network often exhibits multiple stable points, bifurcations, and oscillations. The most important implication of noise in critical cellular processes is that, in spite of identical initial conditions, with time, different cells may evolve along distinct pathways.33 In this context, the presence of noise raises a different set of questions concerned with the robustness, reliability, and sensitivity of the P450 cycle. For example, how does the cycle accurately respond to the environmental changes (like that of the temperature or concentrations) even though some intermediate molecules behave in a probabilistic manner? In order to investigate whether the stochastic effects or noise will drive the system randomly into bifurcation or a certain branching pathway, another modeling attempt concerned with the branching reactions was also

performed. Considering that H2O and H2O2 are products of the two branching ways in the cycle, their simulation results are displayed to analyze and illustrate the behavior of the branching ways of the system. It is noteworthy that the P450 cycle is very complex and not fully understood. Defining the model has involved making simplified assumptions that reduce complex processes into simpler ones that can still represent the biological processes well enough to explain the observed data. When considering the P450 cycle, we notice that one branching point exists in this system. Therefore, the model proposed in Scheme 1 may be somewhat oversimplified for the branching reactions. In order to investigate the reactions more clearly, in the following model (reactions with an asterisk in Table 1), the stable complex intermediate of Fe3+-OOH- 28 was introduced between ES4 and ES5 in Scheme 1 (as EC shown in Figure 7). The reactions concerned with this intermediate are the following: k10

Fe3+-O22-‚RH(ES4) + H + 98 Fe3+-OOH-‚RH(EC) k11

Fe3+-OOH-‚RH(EC) + H+ 98 Fe-O3+‚RH(ES5) + H2O k12

Fe3+-OOH-‚RH(EC) + H + 98 Fe3+‚RH(ES1) + H2O2 In this modeling attempt, all three second-order reactions were set at 6000 µM-1 min-1, a reasonable rate for diffusioncontrolled encounters between large and small molecules, with pH 7 ([H+] ) 100 nM), and the other conditions were the same as those in Figure 1. Figure 7 depicts the modeling results and reveals that, in essence, there is no significant difference in the whole trend of the system as compared to Figure 1, and neither is even for the H2O and H2O2. Meanwhile, no periodic oscillations were found in this model. These observations

Cytochrome P450 Catalytic Cycle

Figure 8. Reaction loops analyzed in the stochastic modeling. Loop 1 represents the series of reactions producing ES1 from EC by way of EC f ES1 f ES2 f ES3 f ES4 (f EC), while loop 2 represents the series of reactions producing ES5 from EC by way of E f ES1 f ... f EC f ES5 f ... (f E).

suggest that the original model of Scheme 1 for analyzing the P450 cycle is, though simplified, still reasonable. In order to further identify whether a dependence of the bifurcation on the individual rate constants (k11, k12) exists, k11 and k12 constants were then chosen to vary to compute the time evolution of H2O2 and H2O and the other conditions were kept the same as those in Figure 1. By trying different combinations of the k11 and k12 values, we found that, when k11 ) 6000 µM-1 min-1, if k12 e 200 µM-1 min-1, almost only loop 2 in Figure 7 involving the production of H2O of the system persists, whereas the formation of H2O2 can seldom be observed. Likewise, for k12 ) 6000 µM-1 min-1, if k11 e 500 µM-1 min-1, it is the loop producing the H2O2 that mainly presents (loop 1 in Figure 8) and loop 2 can hardly proceed. These results were based on the average results of 10 runs of the respective conditions. The observations show that for a stochastic simulation it may in principle be possible for the populations to be extinct for a greater propensity of reactions, since a molecule reacts with its counterpart with the same probability, irrespectively of how they are located in space. However, in a deterministic simulation like ODEs, without specifying k12 ) 0 or k11 ) 0 µM-1 min-1, theoretically, it is impossible to cease their corresponding branch reactions. If this modeling is applicable to the P450s, such occasions may occur that changes in each of the two steps will produce crucial impacts on the dynamic behaviors of the system, which may be an interesting indication for further experimental study. From the above results, an interesting fact is discovered that the two limits of k11 and k12 are different, that is, k11 ) 500 µM-1 min-1, whereas k12 ) 200 µM-1 min-1 (when their counterparts remain 6000 µM-1 min-1). As we all know, typically, a stochastic simulation of intracellular reactions is carried out in a model system whose reactants are uniformly distributed in space. This means that a molecule is assumed to have the same probability to interact with any of its reaction parents in a stochastic modeling framework. Therefore, seemingly, the limit values of k11 and k12 should be the same, since they are both the branching reactions of k10 and seem to share

J. Phys. Chem. B, Vol. 111, No. 16, 2007 4259 the same chemical environment. However, why are the two limits actually different? A physical interpretation to this phenomenon may be as follows: Loop 2 is a consumption reaction cycle, where the reactant molecules are constantly decreased while the product (P) molecules increased in the cycle. Therefore, only by increasing the propensity of the reactants can the chances for the reaction proceeding be increased when its reaction rate is small (i.e., k11 is less than 500 µM-1 min-1). However, no molecules are consumed in loop 1. Thus, the number of the reactants for this loop will not be decreased. Accordingly, the reaction chances of this loop can still remain large even when its reaction rate is relatively small (when k12 is slightly larger than 200 µM-1 min-1). These results led us to believe that the rates of the two branching reactions, under physiological conditions, are different. The above tries have led to different predictions as to the bifurcations of the system, which may provide guidance for further experimental studies. However, in order to more extensively investigate the general mechanism and behavior from the proposed model, a supplement of more accurate in ViVo or in Vitro experimental values for each reaction parameter is still needed. 4. Conclusions A stochastic version of the cytochrome P450 catalytic cycle was built in the present work, which provides valuable insights into the functioning of the biological system considering stochastic effects. During our modeling practice, strong fluctuations are discovered appearing in the entry along with weak fluctuations presenting in the middle sections of the cycle, which phenomenon, we deduce, is important for the P450 natural performance. In addition, on the basis of a microscopic view of points obtained from the work, we conclude that, basically, the fundamental feature of the P450 cycle is deterministic in nature. The system evolves along a relatively fixed pattern from its initial state, which is possibly related to the robustness property of the P450 enzyme. Meanwhile, the simulation results of the two branching reactions, that is, the formations of H2O and H2O2, indicate that the rates of the two reactions are physiologically different, and changes in the two rates can cause different impacts on the behaviors of the enzyme. These findings will advance the understanding of the catalytic mechanism of cytochrome P450s. Acknowledgment. The authors are grateful to Prof. F. Peter Guengerich for helpful discussions. This work was supported by the Youth Teacher Fund of Dalian University of Technology, the Doctoral Fund of Dalian University of Technology, and the Doctoral Fund of Dalian Fisheries University. References and Notes (1) Guengerich, F. P. Chem. Res. Toxicol. 2001, 14, 611-650. (2) Shaik, S.; Kumar, D.; de Visser, S .P.; Altun, A.; Thiel, W. Chem. ReV. 2005, 105, 2279-2328. (3) Ogliaro, F.; de Visser, S. P.; Cohen, S.; Sharma, P. K.; Shaik, S. J. Am. Chem. Soc. 2002, 124, 2806-2816. (4) Gigon, P. L.; Gram, T. E.; Gillette, J. R. Mol. Pharmacol. 1969, 5, 109-122. (5) Diehl, H.; Scha¨delin, J.; Ullrich, V. Hoppe-Seyler’s Z. Physiol. Chem. 1970, 351, 1359-1371. (6) Gander, J. E.; Mannering, G. J. Kinetics of hepatic cytochrome P-450-dependent mono-oxygenase systems. In Hepatic Cytochrome P-450 Monooxygenase System; Schenkman, J. B., Kupfer, D., Eds.; Pergamon Press: New York, 1982; pp 667-697. (7) Bell-Parikh, L. C.; Guengerich, F. P. J. Biol. Chem. 1999, 274, 23833-23840. (8) Kozuch, S.; Shaik, S. J. Am. Chem. Soc. 2006, 128, 3355-3365.

4260 J. Phys. Chem. B, Vol. 111, No. 16, 2007 (9) Yun, C.-H.; Miller, G. P.; Guengerich, F. P. Biochemistry 2001, 40, 4521-4530. (10) Guengerich, F. P.; Krauser, J. A.; Johnson, W. W. Biochemistry 2004, 43, 10775-10788. (11) Wang, Y.; Li, Y.; Li, Y.; Ma, X.; Yang, S.; Yang, L. J. Phys. Chem. B 2006, 110, 10139-10143. (12) Elowitz, M. B.; Levine, A.; Siggia, E.; Swain, P. Science 2002, 297, 1183-1186. (13) Levin, J. E.; Miller, J. P. Nature 1996, 380, 165-168. (14) Becskei, A.; Serrano, L. Nature 2000, 405, 590-593. (15) Gardner, T.; Cantor,C.; Collins, J. Nature 2000, 403, 339-342. (16) Smith, G. D. Modeling the Stochastic Gating of Ion Channels. In Computational Cell Biology; Fall, C. P., Marland, E. S., Wagner, J. M., Tyson, J. J., Eds.; Springer-Verlag: New York, 2002; pp 285-319. (17) Gillespie, D. T. J. Phys. Chem. 1977, 81, 2340-2361. (18) Gibson, M. A.; Bruck, J. J. Phys. Chem. A 2000, 104, 1876-1889. (19) Kastner, J.; Solomon, J.; Fraser, S. DeV. Biol. 2002, 246, 122-131. (20) Kierzek, A. M.; Zaim, J.; Zielenkiewicz, P. J. Biol. Chem. 2001, 276, 8165-8172. (21) Guengerich, F. P.; Krauser, J. A.; Johnson, W. W. Biochemistry 2004, 43, 10775-10788. (22) Houston, J. B. Biochem. Pharmcacol. 1994, 47, 1469-1479. Radius cell ) 50 uM, assuming 50 pmol of P450 1A2/mg of microsomal protein.

Wang et al. (23) Korzekwa, K. R.; Krishnamachary, N.; Shou, M.; Ogai, A.; Parise, R. A.; Rettie, A. E.; Gonzalez, F. J.; Tracy, T. S. Biochemistry 1998, 37, 4137-4147. (24) Backes, W. L.; Kelley, R. W. Pharmacol. Ther. 2003, 98, 221233. (25) Gillespie D. T. MarkoV Processes: An Introduction for Physical Scientists; Academic Press: San Diego, CA, 1992. (26) Shimada, T.; Yamazaki, H.; Mimura, M.; Inui, Y.; Guengerich, F. P. J. Pharmacol. Exp. Ther. 1994, 270, 414-423. (27) Ortiz de Montellano, P. R.; De Voss, J. J. In Cytochrome P450: Structure, Mechanism, and Biochemistry, 3rd ed.; Ortiz de Montellano, P. R., Ed.; Kluwer Academic/Plenum Publishers: New York, 2005; pp 183245. (28) Davydov, R.; Makris, T. M.; Kofman, V.; Werst, D. E.; Sligar, S. G.; Hoffman, B. M. J. Am. Chem. Soc. 2001, 123, 1403-1415. (29) Griffin, B. W.; Peterson, J. A. Biochemistry 1972, 11, 4740-4746. (30) Mueller, E. J.; Loida, P. J.; Sligar, S. G. In Cytochrome P450: Structure, Mechanism, and Biochemistry, 2nd ed.; Ortiz de Montellano, P. R., Ed.; Plenum Press: New York, 1995; pp 83-124. (31) Guengerich, F. P. J. Biol. Chem. 1991, 266, 10019-10022. (32) Guengerich, F. P.; Johnson, W. W. Biochemistry 1997, 36, 1474114750. (33) Arkin, A.; Ross, J.; McAdams, H. H. Genetics 1998, 149, 16331648. (34) Oprian, D. D.; Gorsky, L. D.; Coon, M. J. J. Biol. Chem. 1983, 258, 8684-8691.