Beyond homogeneity: Assessing the validity of the Michaelis–Menten rate law in spatially heterogeneous environments

The Michaelis–Menten (MM) rate law has been a fundamental tool in describing enzyme-catalyzed reactions for over a century. When substrates and enzymes are homogeneously distributed, the validity of the MM rate law can be easily assessed based on relative concentrations: the substrate is in large excess over the enzyme-substrate complex. However, the applicability of this conventional criterion remains unclear when species exhibit spatial heterogeneity, a prevailing scenario in biological systems. Here, we explore the MM rate law’s applicability under spatial heterogeneity by using partial differential equations. In this study, molecules diffuse very slowly, allowing them to locally reach quasi-steady states. We find that the conventional criterion for the validity of the MM rate law cannot be readily extended to heterogeneous environments solely through spatial averages of molecular concentrations. That is, even when the conventional criterion for the spatial averages is satisfied, the MM rate law fails to capture the enzyme catalytic rate under spatial heterogeneity. In contrast, a slightly modified form of the MM rate law, based on the total quasi-steady state approximation (tQSSA), is accurate. Specifically, the tQSSA-based modified form, but not the original MM rate law, accurately predicts the drug clearance via cytochrome P450 enzymes and the ultrasensitive phosphorylation in heterogeneous environments. Our findings shed light on how to simplify spatiotemporal models for enzyme-catalyzed reactions in the right context, ensuring accurate conclusions and avoiding misinterpretations in in silico simulations.


Introduction
Enzymes play significant roles in regulating the rates of reactions in living organisms, influencing essential processes such as metabolism, signal transduction, and cellular regulation [1][2][3][4].For over a century, the Michaelis-Menten (MM) equation has been the prevailing framework for describing rates of enzyme-catalyzed reactions [5][6][7][8][9][10][11]. Specifically, the MM rate law effectively describes the rate of product (P) accumulation in a single enzyme-catalyzed mechanism in terms of substrate concentration (S): where k cat is the catalytic constant, and K M is the MM constant, and E T is the total enzyme concentration.The MM rate law was proposed by Henri [5] and Michaelis and Menten [6], and derived by Briggs and Haldane using an approach known as the quasi-steady state approximation (QSSA) [7].This approximation has been referred to as the standard quasi-steady state approximation (sQSSA) [11][12][13][14].The sQSSA has also been utilized to describe a variety of observable biomolecular interactions between genes and transcription factors [15,16], and receptors and ligands [17,18].The MM rate law based on the sQSSA is accurate when the enzyme concentration is low enough to be E T � S T + K M , where S T is total substrate and product concentration [11][12][13][14], referred to as the validity condition of the sQSSA throughout this study.Since the concentrations of enzyme and substrate are roughly comparable in protein interaction reactions, the sQSSA model can be inaccurate [11,19].To resolve this inaccuracy of the sQSSA model, an alternative model based on the total quasi-steady state approximation (tQSSA) was proposed [11,13,14,[20][21][22][23].Even when the concentration of enzymes is high, and thus the validity condition of the sQSSA is not satisfied, the tQSSA model is accurate [11,13,14,[20][21][22][23].Thus, the tQSSA model has been recognized as a more reliable approximation tool for the enzyme-catalyzed reaction model than the sQSSA model.
Both the sQSSA and the tQSSA models are based on an ordinary differential equation (ODE) that assumes homogeneous distributions of species, including enzymes and substrates.However, this assumption appears not to be valid within the cell, unlike in in vitro experimental environments.Specifically, enzymes are localized in particular organelles depending on their type and function.For instance, enzymes managing the respiratory chain and oxygen reduction reactions are highly concentrated in mitochondria [24].Furthermore, cytochrome P450 (CYP) enzymes, essential for drug metabolism, are located in microsomes in the endoplasmic reticulum (ER) of liver cells [25].Thus, it would be inaccurate to use an ODE to describe enzyme-catalyzed reactions in cells where the assumption of the homogeneous distribution is violated.Instead, employing a partial differential equation (PDE) that accounts for the spatial distribution of chemical species would be a more valid approach.
Previous studies have shown that the applicability of the sQSSA, and thus the MM rate law, can be extended to PDE when the local concentration can be assessed across the entire spatial domain.Specifically, when E T � S T + K M (the validity condition of the sQSSA in ODE) is satisfied at each point within the domain, using the MM rate law in PDE is accurate under various diffusion conditions: (1) when the diffusion time scale is shorter than the fast reaction time scale or longer than the slow reaction time scale [26], (2) when the diffusion time scale aligns with the slow kinetic time scale [27], and (3) when the dissociation of the complex is significantly faster than the diffusion [28].However, verifying this validity condition necessitates measuring enzyme and substrate concentrations at every point within the entire cell, which is impossible.On the other hand, the spatial average concentrations of enzyme and substrate can be measured using techniques such as western blot [29], UV-Vis spectroscopy [30] and Enzyme-Linked Immunosorbent Assay [31].This necessitates a validity condition based on spatial average concentrations rather than the local concentration.
Here we investigated whether the sQSSA in PDE is valid when where � E T and � S T are the spatial average of the total enzyme concentration and the total substrate and product concentration across the domain, respectively.Specifically, we showed that the sQSSA in PDE can still be accurate although the spatial averages of molecular concentrations do not satisfy � E T � � S T þ K M , when the substrate and the enzyme are localized in different regions and their diffusion is slower than slow reactions.Conversely, we found that the sQSSA in PDE may introduce substantial errors, even if the spatial average concentrations satisfy especially when the substrate and the enzyme are localized in the same region under conditions of slow diffusion.As a result, employing the sQSSA model to simulate drug metabolism in the liver, where the CYP enzyme is localized in the ER, leads to considerable error in predicting drug clearance, the rate of a drug's breakdown in the liver.Unlike the sQSSA, using the tQSSA in PDE was accurate for all of these cases because the tQSSA in ODE is generally valid regardless of enzyme concentration.These findings imply that in heterogeneous spaces, where validating the sQSSA at every point is challenging, utilizing the sQSSA in PDE poses a risk.Thus, the tQSSA should be used to model enzyme-catalyzed reactions in heterogeneous spaces.

QSSAs for ODE describing enzyme-catalyzed reactions
An enzyme-catalyzed reaction can be described by the following ODE system based on massaction kinetics (Fig 1 ): where E and C are the concentration of the enzyme and the complex, respectively.The enzyme reversibly binds the substrate to form the complex and then the complex irreversibly dissociates into product and enzyme.Note that E T � E + C and S T � S + C + P are always conserved.By utilizing conservation laws, this model can be effectively expressed as an ODE system with two variables, S and C, as follows: where P can be simply obtained by P = S T − C − S. We can further simplify the model by assuming that C reaches a quasi-steady state (QSS) more rapidly than S and then both C and S slowly track along the QSS.The approximation to QSS (i.e., QSSA) can be obtained by solving dC/dt = 0: where is the MM constant.By substituting Eq 4 for Eq 3, the MM rate law can be derived: As the MM rate law is derived by using the sQSSA for ODE systems, we refer it to as the sQSSA o model throughout the study.The sQSSA o model is accurate when E T is sufficiently small (i.e., E T � S T + K M ) [11][12][13][14].Thus, when E T and S T are comparable, such as in protein interaction reactions, using the sQSSA o model can lead to significant errors [11,19].This inaccuracy arises from treating S as a slow variable, even though S is also affected by fast binding and unbinding reactions [11,32].In other words, the time scale of S is not significantly longer than that of C, violating the underlying assumption of the sQSSA.Such inaccuracy can be resolved by introducing the total substrate concentration Ŝ ¼ S þ C. By using Ŝ, Eq 3 can be rewritten as follows: As Ŝ is now solely affected by a slow reaction (i.e., d ), the time scale of Ŝ is significantly longer than that of C. Thus, C is more likely to reach QSS before Ŝ changes appreciably.The QSSA of C given Ŝ can be obtained by solving dC dt ¼ 0 in Eq 6: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Substituting this QSSA for C (Eq 7) back into the model (Eq 6) results in the reduced model: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where Ŝ ¼ S T À P. As this model is derived with the total quasi-steady state approximation for ODE systems [11,13,14,[20][21][22][23], we referred it to as the tQSSA o model throughout the study.It was shown that the tQSSA o model is accurate when the condition ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi is always satisfied, it has been believed that the tQSSA o model is generally accurate [11,13,14,[20][21][22][23].In particular, when the concentration of enzymes is high and thus the sQSSA o model is not valid, the tQSSA o model is accurate [11,13,14,[20][21][22][23].Although Eilertsen and Schnell recently identified conditions where even the tQSSA o model becomes inaccurate (i.e., the region outside of [33], this limitation of the tQSSA o is confined to a very small region of the parameter space.In particular, even when S T /K M � 1 and thus the above conditions are not satisfied [33], it was also shown that the tQSSA o model still performs reasonably well.Thus, the tQSSA o model has been recognized as a reliable approximation tool for the enzyme-catalyzed reaction model [11].

QSSAs for PDE describing enzyme-catalyzed reactions
So far, we have discussed the application of QSSAs in ODE systems, assuming the even distribution of species in space.To incorporate the spatial distribution of species, we can employ the following PDE system [26][27][28] defined on a bounded region with a smooth boundary O: where diffusion coefficients of S, E, C, and P are denoted as D S , D E , D C , and D P , respectively.
In this study, we mainly used a one-dimensional domain O = (0, L), where L is the length of the domain.Besides, we used zero-Neumann boundary conditions and continuous initial conditions (ICs) S(0, x) = S 0 (x), E(0, x) = E 0 (x), C(0, x) = C 0 (x), and P(0, x) = P 0 (x).Analogous to the ODE system, we can express Eq 9 as follows by using E T = E + C: Note that E T is no longer pointwise conserved because the diffusion can change the concentration of the enzyme.The above PDE system (Eq 10) can be simplified by applying QSSA as in the ODE system when the time scales are separated.Since species not only react but also diffuse throughout the cell, it is essential to compare the diffusion time scale with the reaction time scale.When the diffusion time scale is comparable to or shorter than the fast reaction time scale, the concentrations of each species are rapidly homogenized, resulting in similar dynamics to ODE at each point [26].Thus, we focused on the situation when the diffusion time scale is comparable to or longer than the slow reaction time scale.Under this condition, C reaches QSS before the concentrations of other species change significantly by diffusion and reaction [27].Then, by using the sQSSA o (CðSÞ ¼ E T S SþK M ), the PDE model (Eq 10) can be reduced as follows [27]: where CðSÞ ¼ E T S SþK M .As the sQSSA o is used in the PDE system, we refer to this model (Eq 11) as the sQSSA p .
As in the ODE (Eq 6), we introduce a slow variable Ŝ ¼ S þ C to the PDE system (Eq 10) as follows:

PLOS COMPUTATIONAL BIOLOGY
Then by using the tQSSA o of Cð ŜÞ (Eq 7), we obtained the following reduced PDE system: where ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi . We refer to this model (Eq 13) as the tQSSA p .

The accuracy of the sQSSA p , but not the tQSSA p varies depending on whether the environment is homogeneous or heterogeneous
We explored the accuracy of the sQSSA p (Eq 11) and the tQSSA p (Eq 13) in both homogeneous and heterogeneous environments, by comparing their behaviors to those of the full model (Eq 9).To simulate these models, we first chose biologically realistic values for the parameters (D S , D E , D C , D P , k f , k b , k cat , and L).For instance, we assigned a diffusion coefficient of D � = 0.2μm 2 /s to all species, which is lower than typical protein diffusion coefficients [34] but corresponds to the diffusion coefficient of large-sized proteins (e.g., PER2 protein complex) [35].Furthermore, since the protein-protein binding occurs with a rate on the order of 10 6 M −1 s −1 , often ranging up to *10 8 M −1 s −1 [36,37], we have used k f = 3.4 � 10 6 M −1 s −1 .We set k b and k cat to 60s −1 and 3.2s −1 , respectively, ensuring that K M , k cat , and k cat /K M are in the range of typical enzyme kinetic parameters [38].Finally, the length of the domain O, L, was set to 30μm, which falls within the human cell size range [39].We first investigated the homogeneous ICs so that the validity condition of the sQSSA o (E T � S T + K M ) is satisfied at every point of the domain ( . Specifically, we increased the concentrations of S and E near x = 30μm and decreased the concentrations in other regions.As a result, the validity condition of the sQSSA o (S T + K M � E T ) is not satisfied near x = 30μm (Fig 2D).In this region, the sQSSA p overestimates P of the full model (Fig 2E).Since the majority of S is localized at this invalid region, the majority of P is generated in that region, and thus, the sQSSA p also overestimates � P compared to the full model (Fig 2F).On the other hand, even with these heterogeneous ICs, the tQSSA p accurately simulates P (Fig 2E and 2F).
The impact of a heterogeneous environment on the accuracy of the sQSSA p may diminish with faster diffusion because the initial spatial distribution of species is homogenized and the spatial averages of the concentrations satisfy the validity of the sQSSA o (Fig 2A).Indeed, as the diffusion coefficient (D � ) increases, the accuracy of the sQSSA p becomes accurate (S1 Fig) .Specifically, when the diffusion time scale (L 2 /D � ) is comparable to the time scale of the slowest reaction (1/k cat ), the sQSSA p accurately captures the behavior of the full model.On the other hand, when the diffusion time scale is *100-fold slower than the time scale of the slowest reaction, the sQSSA p is inaccurate.Taken together, the inaccuracy of the sQSSA p arises when the diffusion is slower than all of the reactions.

PLOS COMPUTATIONAL BIOLOGY
We next investigated the opposite case: can the sQSSA p be accurate, even if � E T � � S T þ K M is not satisfied?For this, we first used the homogeneous ICs where the validity condition of the sQSSA o is not satisfied at any point (Fig 2G).Obviously, the sQSSA p overestimates both P and � P of the full model (Fig 2H and 2I).In contrast, the tQSSA p accurately captures the dynamics of the full model (Fig 2H and 2I).Now, we heterogenized the ICs while maintaining the spatial average concentrations of Fig 2G (Fig 2J).Specifically, by altering the distribution of S and E, we made the validity condition of the sQSSA o be satisfied for x < 15μm, but not for x > 15μm.With these heterogeneous ICs, unexpectedly, the sQSSA p accurately captures the dynamics of the full model (Fig 2K and 2L).This is because most of S is localized where the validity condition of the sQSSA o is satisfied, and thus the majority of P is generated in the valid region (Fig 2J, blue region).As a result, the sQSSA p can accurately simulate the production rate of P. The tQSSA p also accurately captures the dynamics of the full model (Fig 2K and 2L).
Taken together, the validity of the sQSSA p critically depends on the spatial distribution of the ICs.The sQSSA p can be invalid even when the spatial average concentration satisfies the validity condition of the sQSSA o ( � , and conversely, the sQSSA p can be valid even when Therefore, using the validity condition of the sQSSA o based on spatial average concentrations of species is not enough to determine the validity of the sQSSA p .On the other hand, the tQSSA p is accurate regardless of the spatial distribution of the ICs.
Unlike the tQSSA p , the sQSSA p may not be applicable in the in vivo environment where the enzyme is localized In the previous section, we demonstrated that the validity of the sQSSA p cannot be ensured solely by the spatial average concentrations when species are not homogeneously distributed within a domain.This situation is common within cells, as specific enzymes are localized within distinct subcellular organelles such as the nucleus, the ER, lysosomes, and mitochondria (Fig 3A).For instance, in hepatocytes, cytochrome P450 enzymes, which metabolize drugs, are localized in the ER [25].
To mimic in vivo environments with the localized enzyme distribution, we utilized the ICs where the enzyme is localized in a narrow region near x = 15μm (Fig 3B, dashed line), while the substrate is uniformly distributed throughout the entire cell (Fig 3B , solid line).The ICs do not satisfy the validity condition of the sQSSA o in a narrow region near x = 15μm.As a result, the sQSSA p model overestimates the production of P (Fig 3C and 3D).This is because the majority of the P is produced in the region where the enzyme is localized, and thus sQSSA o is not valid.However, the tQSSA p model accurately captures the dynamics of the full model (Fig 3C and 3D).Taken together, it is safer to use the tQSSA p model than the sQSSA p model in an in vivo environment.
In vitro experiments are usually performed with enzyme and substrate concentrations similar to those in in vivo experiments.To mimic this in vitro environment, we homogenized the ICs in Fig 3B while maintaining the spatial average concentrations (Fig 3E).With this homogenization, the validity condition of the sQSSA o (E T � S T + K M ) is satisfied throughout the domain.As a result, not surprisingly, in this case, both the sQSSA p and tQSSA p models accurately approximate the dynamics of the full model (Fig 3F and 3G).This indicates that although the sQSSA p model provides a precise approximation of the full model in in vitro, it cannot be directly translated to the in vivo environment.This is consistent with previous studies showing that utilizing sQSSA o to extrapolate in vitro drug clearance to in vivo drug clearance results in a substantial overestimation of the drug clearance, while the tQSSA o model provides a reliable estimate for drug clearance [40,41].Furthermore, in in vivo environments, directly measuring enzyme concentrations throughout cells (e.g., CYP in hepatocytes) is challenging.In the absence of such information, one can investigate how drug clearance changes depending on the heterogeneity of the enzyme distribution in the cell.As the heterogeneity of the initial distribution of E increases (Fig 3H), the initial velocity of production formation considerably decreases (Fig 3I).This is accurately captured by the tQSSA p model, but not the sQSSA p model (Fig 3I).This indicates that the sQSSA p model tends to overestimate drug clearance to a larger extent as the cellular environment becomes more heterogeneous.

Unlike the tQSSA p , the sQSSA p can create artificial ultrasensitivity and spatial patterns
So far, we have investigated how the accuracy of QSSAs can vary in a simple enzyme-catalyzed model.Here, we examined biological systems that encompass multiple enzymes, which can exhibit complex behaviors, such as ultrasensitivity, bistability, and oscillation [42].One prominent example of such reactions is the Goldbeter-Koshland (GK) mechanism (Fig 4A).The GK mechanism describes two interlocked enzyme-catalyzed reactions with a single substrate [11,43].Specifically, the substrate (S) is phosphorylated by the kinase (E), and the phosphorylated substrate (S P ) is dephosphorylated by the phosphatase (D).We extended the GK mechanism to the PDE system to represent the spatial distribution of species within the cell in a two-dimensional domain (O = (0, L) × (0, L)) and derived its sQSSA p and tQSSA p models (Fig 4A ) (see Methods for details).
We explored whether the sQSSA p and the tQSSA p models capture the behaviors of the full model across varying concentrations of substrate and enzymes.Specifically, we divided the domain into two regions and gave them two different ICs: one with a high S T /D T ratio (y � 15μm), and another with a small S T /D T ratio (y < 15μm) (Fig 4B Ultrasensitivity is important for pattern formation across various biological and chemical reactions such as cell division [42,[44][45][46].We examined how this 'artificial' ultrasensitivity caused by the sQSSA p affects pattern formation.For this, we used ICs that do not satisfy the validity condition of the sQSSA o across all regions (i.e.,  around one cause significant changes in Ŝp =S T due to this artificial ultrasensitivity (Fig 4C (iii)).Consequently, a vertical striped pattern emerges and combines with the horizontal pattern, resulting in the observed grid-like pattern (Fig 4E (iii)).This indicates that using the MM equation can generate misleading artificial ultrasensitivity and patterns, although it has been frequently utilized to depict pattern formation [47,48].Thus, utilizing the tQSSA is a safer option to investigate pattern formation.

Discussion
Although the MM rate law has been the prevalent means for describing enzyme-catalyzed reactions for over a century [5][6][7][8][9][10][11], the applicability of the MM rate law remains elusive when species are heterogeneously distributed.However, our study demonstrated that the validity condition of the sQSSA p (i.e., the MM rate law in a spatiotemporal system) is different from that of the sQSSA o (Fig 2).Specifically, even if the spatial average concentrations of species satisfy the validity condition of the sQSSA o , the sQSSA p model can be inaccurate.And conversely, the sQSSA p model can be accurate even if the spatial average concentrations do not satisfy the validity condition of the sQSSA o .Thus, the validity condition based on spatial average concentration does not ensure the accuracy of the sQSSA p model.Therefore, to assure its accuracy, it is necessary to verify that ICs satisfy the validity condition of the sQSSA o at each point in the domain.However, this is not feasible because it is challenging to measure concentrations at every point in the cell.Thus, since the tQSSA p is valid in most cases, our findings suggest that the tQSSA p offers the best way to reduce and investigate spatiotemporal models.
For systems with fast diffusion, the dynamics of PDE and ODE models are very similar because the concentrations of each species homogenize quickly [26].In such systems, the validity condition of the sQSSA p would be similar to that of the sQSSA o .However, we found that for systems with slow diffusion, the validity of the sQSSA p becomes different from the sQSSA o (Figs 2 and 4).This indicates that in models involving slow diffusion, such as 1) measuring dopamine concentration in the brain [49] and 2) detecting current differences with biosensors [50][51][52], the MM rate law could be inaccurate even when the validity condition of the sQSSA o is satisfied.Thus, when reducing spatiotemporal models with slow diffusion, it is safer to use the tQSSA p than the sQSSA o to avoid potential inaccuracies.
Drugs administered into the body are primarily broken down by CYP enzymes located in liver cells [25].Drug clearance serves as a key pharmacokinetic parameter for predicting human systemic drug disposition, as cited in over 65,000 papers [40,42,[53][54][55].While the MM rate law has conventionally been used to estimate drug clearance, its validity in drug clearance has only recently been verified [40,41].Recent studies suggested that this approach is only valid when the spatial average CYP concentration is significantly smaller than the MM constant of the drug (i.e., � E T =K M < 0:1) [40,41].Otherwise, the tQSSA model should be used to predict the drug clearance [40,41].However, our findings indicate that this criterion for using the MM rate law, based on the spatial average concentration, can be misleading because the liver cell is heterogeneous (Fig 3).Therefore, even when � E T =K M < 0:1 is satisfied, utilizing the tQSSA p may be a safer approach for predicting drug clearance.
In experimental settings including in vivo environments, accurately measuring concentrations of substances within cells poses a challenge, which hampers the investigation of how spatial differences affect enzyme-catalyzed reactions.While utilizing various initial conditions can represent different levels of heterogeneity, simulating the full model, encompassing all binding and unbinding interactions, is computationally expensive.Therefore, employing simplified models using sQSSA p or tQSSA p is more efficient.However, the sQSSA p model may not accurately reflect the behavior of the full model in the presence of spatial heterogeneity.For example, our simulation showed that the sQSSA p model underestimated the impact of spatial heterogeneity on the initial velocity of enzyme-catalyzed reactions (Fig 3H and 3I).Thus, using the tQSSA p model, which accurately captures the behavior of the full model regardless of spatial heterogeneity, is a reasonable solution.Indeed, the tQSSA p model accurately captures the range of the initial velocities of enzyme-catalyzed reactions depending on the spatial heterogeneity of enzyme distributions (Fig 3I).
We have investigated whether QSSAs in the PDE system provide accurate approximations when species are heterogeneously distributed (Figs 2-4).Our work shows that the accuracy of the tQSSA p model, but not the sQSSA p model, is not affected by the heterogeneity level.However, to determine whether this result is valid universally, it would be essential to investigate this with a broader range of conditions.Furthermore, our findings are based on numerical simulations rather than theoretical derivations.Notably, recent theoretical studies have identified a small region where the tQSSA o becomes inaccurate, although the tQSSA o has been believed to be generally accurate [33,56].Therefore, the theoretical investigation of the validity of the tQSSA p should be a crucial focus for future research.Furthermore, a recent study proposed a more accurate model than the MM rate law by using the effective time-delay scheme [57].Investigating this new scheme in heterogeneous environments would be also interesting in future work.
Even when the tQSSA p accurately approximates the original PDE systems, the tQSSA p may not be accurate in a stochastic context.Using the tQSSA o for stochastic simulations is usually accurate [11,[58][59][60][61][62][63][64] and thus has been widely used [65][66][67][68].However, a recent study has shown that using the tQSSA o for stochastic simulations can be inaccurate even though it is accurate in a deterministic simulation in which the molar ratio of E and S is close to 1:1 and their binding is tight [69].This implies that the validity condition of the tQSSA o in stochastic systems is different from that in deterministic systems.Therefore, it would be interesting in future work to investigate the validity of the tQSSA p in deterministic and stochastic contexts.

Numerical simulation for a simple enzyme-catalyzed reaction model (1-D model)
In Figs 2 and 3, to simulate the enzyme-catalyzed reaction, we numerically solved PDE systems (Eqs 9, 11 and 13) in a one-dimensional space O = (0, L).Let x j = jΔx = jL/(N x − 1) be the N xequispaced grids in O, which are referred to as the Fourier collocation points, to implement the cosine-based spectral method.The time step Δt was defined as Δt = T/N t , where T is the final time and N t is the total number of time steps.We denote the numerical approximation of S(nΔt, x), E(nΔt, x), C(nΔt, x), P(nΔt, x) by S n , E n , C n , P n , respectively.
To facilitate simulation, we split PDE systems into diffusion and reaction parts using the operator splitting method [70].That is, for each time step of the simulation, the diffusion and reaction were calculated separately.We first calculated the diffusion parts in Eqs 9, 11, and 13 as follows:

PLOS COMPUTATIONAL BIOLOGY
After obtaining S * , E * , C * , and P * by calculating the diffusion part, the reaction parts were calculated to obtain S n+1 , E n+1 , C n+1 , and P n+1 .Specifically, for the full model (Eq 9), the following reaction equations were used.
For the sQSSA p model (Eq 11), we defined S * þK M to calculate the reaction parts.Then, we obtained S n+1 and P n+1 as follows: Furthermore, using S n+1 and P n+1 , we obtained ðS nþ1 þK M Þ , and For the tQSSA p model 13), we defined ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi � to calculate the reaction part.Then the following reaction equations were used to obtain Ŝnþ1 and P n+1 .
Then, using Ŝnþ1 and P n+1 , we obtained ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi � � .Then, the heterogeneity of the E 0 was defined as follows: where σ is the standard deviation of the normal distribution, and we set s max ¼ 25mm since E becomes almost constant across the space when σ � 25μm.Additionally, we measured the initial velocity of the enzyme-catalyzed reaction as the average reaction velocity over the first 12 seconds from the initial state, by using:

PDE model describing the GK mechanism
In Fig 4, we extended the GK mechanism [11,43], which consists of two interconnected, single-substrate, enzyme-catalyzed reactions, to the PDE model as follows: where δ S , d S p , δ ES , and d DS p denote the diffusion coefficients of S, S p , ES, and DS p , respectively.
In this model, S T = S + S p + ES + DS p , E T = E + ES, and D T = D + DS p are concentrations of total substrate and product, total kinase, and total phosphatase, respectively.This full model (Eq 14) can also be reduced by applying the sQSSA p , as done in the simple enzyme model (Fig 4A ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi , respectively.Then we obtain the following reduced model of the GK mechanism: Similar to Figs 2 and 3, we utilized the operator splitting method and spectral method.Therefore, the diffusion and reaction were still calculated separately for each time step.We first calculated the diffusion parts in Eqs 14-16 as follows: For the tQSSA p model (Eq 16), we defined ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi . Then, we obtained Ŝnþ1 ; Ŝnþ1 p as follows: For simulations, we used k fe = k fd = 2.22μM The computational codes for implementing all of the models used in this study can be found at https://github.com/Mathbiomed/PDE_QSSA.

Fig 1 .
Fig 1. QSSAs for an enzyme-catalyzed reaction.A model describing an enzyme-catalyzed reaction based on massaction kinetics (Eqs 2 and 9) can be simplified using either the sQSSA o (Eqs 5 and 11) and the sQSSA p (Eq 11) or the tQSSA o (Eq 8) and the tQSSA p (Eq 13).Both approximations of QSS are defined in terms of S or Ŝ^¼ S þ C. https://doi.org/10.1371/journal.pcbi.1012205.g001 Fig 2A).In this case, as expected, both the sQSSA p and the tQSSA p accurately capture P of the full model at each point (Fig 2B), and thus the spatial average of P, � P (Fig 2C).Next, we made the ICs heterogeneous while maintaining the spatial average concentrations of Fig 2A, so that the spatial average concentrations satisfy the validity condition of the sQSSA o

Fig 2 .
Fig 2. The accuracy of the sQSSA p and the tQSSA p in heterogeneous environments.(A) Homogeneous ICs that satisfy the validity condition of the sQSSA o throughout the domain (blue region).E 0 � 2μM, S 0 � 2μM and K M = 18.5μM.(B-C) Both the sQSSA p and the tQSSA p accurately capture the production of P throughout the domain (B) and the spatial average � P � (C).(D) The ICs in (A) were heterogenized so that the validity condition of the sQSSA o model does not hold near x = 30μm (red region), unlike the other region (blue region).Here, E 0 (x) = S 0 (x) = 18.5(tanh(20x/3 − 190) + 1) + 10 −4 μM, which maintain the spatial average concentration with the ICs in (A), was used.(E-F) The tQSSA p model, but not the sQSSA p model, is accurate.(G) The homogeneous ICs that do not satisfy the validity condition of the sQSSA o through the domain.Here, E 0 � 400μM, S 0 � 300μM, and K M = 18.5μM.(H-I) The tQSSA p model, but not the sQSSA p model, is accurate.(J) The ICs in (G) were heterogenized so that the validity condition of the sQSSA o is satisfied in x < 15μm (blue region), but not in x � 15μm (red region).Here, E 0 (x) = 300(tanh(2x/3 − 10) + 1) + 100μM and S 0 (x) = 290 (tanh(10 − 2x/3) + 1) + 10μM, which keep the spatial average concentration with the ICs in (G), were used.(K-L) Both the sQSSA p and tQSSA p models are accurate.For all simulations, C 0 � 0μM, P 0 � 0μM and k f = 3.4 � 10 6 M −1 s −1 , k b = 60s −1 , k cat = 3.2s −1 , D E = D S = D C = D P = D � = 0.2μm 2 /s.https://doi.org/10.1371/journal.pcbi.1012205.g002

Fig 3 .
Fig 3.The sQSSA p , but not the tQSSA p , poses a risk when the enzyme is localized within cells, unlike in an in vitro experiment.(A) Enzymes localized within cellular organelles.(B) The heterogeneous ICs where the sQSSA o are invalid near x = 15μm (red region) but valid in the other region (blue region).Here S 0 � 39 μM, E 0 (x) = 5 � f(x)μM, where f(x) is the normalized probability density function of the normal distribution with the mean of 15 μm and the standard deviation of 0.2 μm.(C-D) Unlike the sQSSA p model, the tQSSA p model accurately captures the production of P throughout the domain (C) and its spatial average � P � (D).(E) The ICs in (B) were homogenized so that the validity condition of the sQSSA o is satisfied throughout the domain while maintaining the spatial average concentration.Specifically, S 0 � 39 μM and E 0 � 5μM were used.(F-G) Both the sQSSA p and tQSSA p models are accurate.(H) Enzyme distributions with varying heterogeneity were constructed using E 0 (x) = 5 � f(x|σ) μM, where f(x|σ) is the normalized probability density function of the normal distribution with the mean of 15μm and the standard deviation of σ.As σ decreases, the heterogeneity increases (see Method for details).(I) For the enzyme distribution with varying heterogeneity, the tQSSA p model, but not the sQSSA p model, accurately captures the initial velocity (I).For all simulations, we used k f = 6.7 � 10 5 M −1 s −1 , k b = 0.53s −1 , k cat = 0.13s −1 , K M = 1μM, and the diffusion coefficients: D E = D C = 0μm 2 /s, and D P = D S = 0.2μm 2 /s.Some parts of Fig 3 were retrieved from Biorender.https://doi.org/10.1371/journal.pcbi.1012205.g003 (i)).As a result, the ICs satisfy the validity condition of the sQSSA o in y � 15μm (i.e., D T K MD þS T � 0:02 � 1), but not in y < 15μm ( D T K MD þS T � 1).In y � 15μm, as E T /D T increases in a horizontal direction (Fig 4B(ii)), the full model exhibits zero-order ultrasensitivity of phosphorylation at the steady state (Fig 4C(i)).That is, at E T / D T < 1 (x < 15μm), the majority of S remains unphosphorylated, and thus the steady-state fraction of total phosphorylated substrate ( Ŝp =S T ) is close to zero (Fig 4C(i)).On the other hand, when E T /D T surpasses 1 at x = 15μm, the phosphorylation of S becomes suddenly predominant, causing Ŝp =S T to sharply increase and approach 1 (Fig 4C(i)).This ultrasensitivity is successfully captured by both the tQSSA p (Fig 4C(ii)) and the sQSSA p models (Fig 4C(iii)).In y < 15μm, the full model does not exhibit ultrasensitivity.That is, Ŝp =S T gradually increases as E T /D T increases and becomes saturated near x = 15μm, where E T /D T � 1 (Fig 4C (i)).Thus, Ŝp =S T does not sensitively change around E T /D T � 1.The tQSSA p accurately captures the smooth changes in phosphorylation of the full model (Fig 4C(ii)).However, the sQSSA p exhibits an 'artificial' ultrasensitivity (Fig 4C(iii)).That is, using the sQSSA p to describe the GK mechanism can result in an 'artificial' ultrasensitivity.
D T K MD þS T � 1) (Fig 4D(i)).Furthermore, we made ICs of S T /D T and E T /D T that have horizontal and vertical striped patterns, respectively (Fig 4D(i) and 4D(ii)).Because E T /D T changes slightly around one along the x-direction, the full model simulates minimal change of Ŝp =S T along the x-direction (Fig 4E(i)).On the other hand, since S T /D T changes considerably along the y-direction, Ŝp =S T changes considerably along the y-direction.As a result, the full model exhibits a horizontally striped pattern (Fig 4E(i)).This is accurately captured by the tQSSA p model (Fig 4E(ii)).However, the sQSSA p model exhibits an artificial grid pattern (Fig 4E(iii)).This artifact arises due to the artificial ultrasensitivity inherent to the sQSSA p .In other words, even slight fluctuations in E T /D T

Fig 4 .
Fig 4. The sQSSA p generates false patterns due to artificial ultrasensitivity.(A) The model diagram depicting the GK mechanism.The full model (Eq 14) based on mass-action kinetics can be reduced by replacing ES and DS p with either the sQSSA p (Eq 15) or the tQSSA p (Eq 16).(B) The heterogeneous ICs where the sQSSA o is valid in y > 15μm (black triangle), but invalid in y < 15μm (white triangle).For this, we used S 0 (x, y) = 500 tanh(3.3y− 50) + 520μM, E 0 (x, y) = 10 tanh(0.2x− 3) + 20μM, and D 0 � 20μM.(C) In y > 15μm, the ultrasensitivity of the full model (i) is accurately captured by both the tQSSA p (ii) and sQSSA p models (iii).On the other hand, in y < 15μm, the sQSSA p model generates artificial ultrasensitivity (iii) unlike the full (i) and tQSSA p models (ii).Here, Ŝ^p =S T ¼ S p =S T for the sQSSA p model, because DS p is assumed to be negligible.(D) Complex ICs where S T /D T exhibits horizontal stripes and E T /D T exhibits vertical stripes.These ICs do not satisfy the validity condition of the sQSSA o throughout the domain.For this, we used S 0 (x, y) = 40 cos(2y/3) + 100 μM, E 0 (x, y) = 5 cos(2x/3) + 100μM, and D 0 � 100 μM.(E) The full model (i) and the tQSSA p model (ii) exhibit horizontally striped patterns.In contrast, the sQSSA p model (iii) results in a grid pattern due to the artificial ultrasensitivity.(C) and (E) were obtained when t = 37.5s.For initial conditions of other species, S p,0 � 0 μM, ES 0 � 0μM, and DS p,0 � 0 μM were used.In addition, we used k fe = k fd = 2.22 � 10 6 M −1 s −1 , k be = k bd = 1.84s −1 , k e = k d = 0.38s −1 , K ME = K MD = 1μM, and d S ¼ d S p ¼ d ES ¼ d DS p ¼ 0:2mm 2 =s.https://doi.org/10.1371/journal.pcbi.1012205.g004
).In particular, ES and DS p are assumed to reach each QSS, resulting in QSSA ofESðSÞ ¼ E T SSþK ME , and DS p ðS p Þ ¼ in the rates of phosphorylation and dephosphorylation in terms of E T , D T , S, and S p : D T S p S p þK MD , where K ME = (k be + k e )/k fe , K MD = (k bd + k d )/k fd .This PLOS COMPUTATIONAL BIOLOGY results ES and Ŝp ¼ S p þ DS p and apply the tQSSA p by substituting ES and DS p with ESð ŜÞ ¼ 1 2 S n Dt ¼ d S DS * ; After obtaining S*, S * p , ES*, and DS * p by calculating the diffusion part, the reaction parts were calculated to obtain S n+1 , S nþ1p , E n+1 , D n+1 , ES n+1 , and DS nþ1 p .Specifically, for the full model (Eq 14), the following reaction equations were used.S nþ1 À S * Dt ¼ À k fe S * E n þ k be ES * þ k d DS p * ; Dt ¼ À k fd S * p D n þ k bd DS * p þ k e ES * ; E nþ1 À E n Dt ¼ À k fe S * E n þ ðk be þ k e ÞES * ; D nþ1 À D n Dt ¼ À k fd S * p D n þ ðk bd þ k d ÞDS * p ; ES nþ1 À ES * Dt ¼ k fe E n S * À k be ES * À k e ES * ; p :In the sQSSA p model (Eq 15), we definedE * T ¼ E n þ ES * , D * T ¼ D n þ DS * p , ES * ðS * Þ ¼ * Dt ¼ À k e ES * ðS * Þ þ k d DS * p ðS * p Þ; p Dt ¼ k e ES * ðS * Þ À k d DS * p ðS * p Þ: