## Abstract

Schlögl's model is the canonical example of a chemical reaction system that exhibits bistability. Because the biological examples of bistability and switching behaviour are increasingly numerous, this paper presents an integrated deterministic, stochastic and thermodynamic analysis of the model. After a brief review of the deterministic and stochastic modelling frameworks, the concepts of chemical and mathematical detailed balances are discussed and non-equilibrium conditions are shown to be necessary for bistability. Thermodynamic quantities such as the flux, chemical potential and entropy production rate are defined and compared across the two models. In the bistable region, the stochastic model exhibits an exchange of the global stability between the two stable states under changes in the pump parameters and volume size. The stochastic entropy production rate shows a sharp transition that mirrors this exchange. A new hybrid model that includes continuous diffusion and discrete jumps is suggested to deal with the multiscale dynamics of the bistable system. Accurate approximations of the exponentially small eigenvalue associated with the time scale of this switching and the full time-dependent solution are calculated using Matlab. A breakdown of previously known asymptotic approximations on small volume scales is observed through comparison with these and Monte Carlo results. Finally, in the appendix section is an illustration of how the diffusion approximation of the chemical master equation can fail to represent correctly the mesoscopically interesting steady-state behaviour of the system.

## 1. Introduction

Molecular reactions are the chemical basis of cellular biological functions. The complex biochemical reaction networks in terms of proteins and macromolecular complexes play a wide range of important roles in cellular signalling, control and regulation (Zhu *et al.* 2004). In recent years, a great amount of attention has been given to the formulation and analysis of models describing reaction networks. We are still looking for new concepts and ideas to illustrate and quantify what is important in these systems.

In this paper, we wish to use the example of a simple model first studied by Schlögl (1972) to illustrate and interlace some of these recent ideas. The major issues are as follows.

### 1.1 Stochastic modelling

Chemical systems inside a cell (especially those of signalling networks involving transcription regulation, protein phosphorylation and GTPases) often involve a small number of molecules of one or more of the reactants. Thus, the traditional method of modelling systems by describing concentration changes with ordinary differential equations (ODEs) and the law of mass action (Zheng & Ross 1991; Murray 2002) is inappropriate. A Markov chain (or master equation) model accounts for the discrete, probabilistic nature of the chemical reactions at the molecular level, but can be difficult to analyse. Diffusion (Fokker–Planck) approximations to the master equation were first developed by Van Kampen (1981) and shown by Kurtz (1976, 1978) to match the solution to the master equation in the thermodynamic limit for some finite time. However, unless the steady state is unique in the macroscopic description (Mansour *et al.* 1981), the two models can disagree in the infinite time limit (see appendix A). Microscopic simulations have validated the master equation as the most accurate description of a reactive process; see Baras & Mansour (1997) for an up-to-date review. In this paper, we illustrate the inability of a deterministic model to capture the behaviour of Schlögl's model. We suggest a new method of analysing the master equation through its spectrum and discuss which model formulation best represents the complete dynamics of the system.

### 1.2 Functional cellular attractors and multiple time scales

In terms of the chemical master equation (CME) formulation (Delbrück 1940; McQuarrie 1968; Arkin *et al.* 1998; Gillespie 2007; Beard & Qian 2008), each stable steady state of the deterministic model corresponds to a peak (while an unstable steady state corresponds to a trough or a saddle) in the stationary probability distribution (the ‘steady-state’ solution to the master equation, which is generally unique Schnakenberg 1976). We refer to these states as functional cellular attractors (FCA) because from the standpoint of cell biology, which is the most relevant field for our theory, they are the states in which the system is most likely to be found and in which the system performs its function(s). In fact, Schlögl's model exhibits multiple time scales: a fast scale where the system relaxes to one of the FCA, and a slow scale over which the system transitions from one FCA to another (Matheson *et al.* 1974). A key question is how long the system remains in each of the FCA and, as we will show, the CME and Fokker–Planck descriptions can yield conflicting answers to this question. Although we restrict ourselves to the case of well-mixed reactions, spatial considerations are a key component of intracellular modelling (Antoine & Lemarchand 2007). See Elf & Ehrenberg (2004) for an analysis of biological bistability including spatial domains.

### 1.3 Thermodynamics of open chemical systems

In the environment of a living cell, biochemical reaction systems operate in an ‘open’ setting (Nicolis & Prigogine 1977; Qian 2006, 2007; Ross 2008). That is, there is a flux of material and/or energy acting on the system. We can no longer rely on the traditional theory of equilibrium thermodynamics to create an accurate model of these systems. The non-equilibrium theory allows the possibility of multiple steady states and non-zero steady-state flux (Kjelstrup & Bedeaux 2008), whereas a closed molecular system tends to a thermal, chemical equilibrium, which is unique and in which the fluxes are zero (Lewis' 1925 principle of detailed balance). There is also a non-zero entropy production rate, which characterizes the non-equilibrium steady state (NESS). Recent developments in the area of fluctuation theorem (Andrieux & Gaspard 2007; Sevick *et al.* 2008) have highlighted the importance of entropy production and its relationship with the irreversible nature of a system. Han & Wang (2007, 2008) have related entropy production (or ‘dissipation cost’) to the robustness of a network. Can entropy production rates predict the relative stability of the FCA? The issues related to this have been controversial in the literature (Callen 1957; Glansdorff *et al.* 1974; Keizer & Fox 1974; Jaynes 1986; Qian 2002; Dewar 2005; Martyusheva & Seleznev 2006), and we seek insights through this concrete example.

This paper is organized as follows. In §2, we introduce the deterministic and stochastic formulations of Schlögl's model. In §3, the equilibrium and NESSs are examined, the region of bistability is discussed and the flux, chemical potential and entropy production rate are defined for both the models. A new way to formulate the model in the case of bistable behaviour is suggested. In §4, a numerical method for the time-dependent and steady-state behaviours is examined, and compared with Monte Carlo simulations and asymptotic results. Section 5 gives a review of the results and an outline of the future work. Appendix A illustrates how the Fokker–Planck approximation to the master equation can fail to accurately predict the correct steady-state distribution.

## 2. The models

We will focus on an autocatalytic, trimolecular reaction scheme, first proposed by Schlögl (1972),(2.1)and(2.2)Let *a* and *b* denote the concentrations of the chemicals A and B, respectively, which are parameters of the system. Let *x* be the concentration of the dynamic chemical X. The system is assumed to be homogeneous in space, and the volume of the system will be denoted as *V*. Owing to the fixed *a* and *b*, the system is not closed, but rather open, with an exchange of chemicals between two material baths. When the chemical potentials (see §3.2) of the two baths are equal, the system reaches chemical equilibrium, as predicted by Gibbs' grand canonical theory. When the baths are not equal, the reaction system is driven.

The deterministic model based on the law of mass action is a first-order, nonlinear ODE (Murray 2002),(2.3)Depending on the parameters, there can be one, two or three non-negative steady states. We are particularly interested in the set of parameters for which there are two stable steady states separated by an unstable steady state, the bistable case.

The stochastic model is given in the form of the CME, an infinite system of coupled ODEs (Beard & Qian 2008). Let integer random variable *n*_{X}(*t*) be the number of X molecules at time *t*, *n*_{A} and *n*_{B} be the number of A and B molecules, which are fixed for a fixed *V*, and *p*_{n}(*t*) be the probability of having *n* X molecules at time *t*: *p*_{n}(*t*)=Pr{*n*_{X}(*t*)=*n*}. The equations are(2.4)and(2.5)where and . The CME may also be referred to as a birth and death process in the theory of Markov processes (Gardiner 1985), where *λ*_{n} and *μ*_{n} are the birth and death rates, respectively.

The rates in the stochastic model are related to the rates in the deterministic model by a factor of *V* which depends on the number of reactants involved in the *i*th reaction. For a reaction involving *m* reactants, the relation will be . In the following analysis, when the volume varies, we consider the ODE reaction rates *k*_{i} and the parameters *a* and *b* to be constant (as opposed to fixing the , *n*_{A} and *n*_{B}). This corresponds to the thought experiment in which the concentrations are fixed while the system volume and numbers of molecules vary correspondingly. Thus, the birth and death rates are written in terms of the ODE parameters as(2.6)and(2.7)

The steady-state probability distribution of the stochastic model can be found by noting that when the system is stationary. We shall call this equation *mathematical detailed balance* following the theory of Markov processes (Schnakenberg 1976; Gardiner 1985; Jiang *et al.* 2004), in which the probability flux in the forward direction is equal to the probability flux in the backward direction at each state in the CME. This leads to the well-known stationary probability,(2.8)Note that the steady state of the stochastic model (i.e. mesoscopic view) is unique for all values of the parameters. The stable and unstable steady states of the ODE are represented as minima and maxima, respectively, in the steady-state probability distribution of the CME.

## 3. Analysis

In this section, we analyse and compare the deterministic (ODE) and the stochastic (CME) models. Under equilibrium conditions, the two models are in complete agreement. However, for the parameters that yield a bistable system, there are stark differences in the predicted behaviour, and in the definition of thermodynamic quantities which describe the system.

### 3.1 Chemical flux and steady-state behaviour

There are four chemical reactions, each of which has a *chemical flux*, which describes how that reaction contributes to the change of *x*:(3.1)The ODE model can be expressed as the sum of the forward fluxes minus the sum of the backward fluxes,(3.2)In the macroscopic setting, the steady state is reached when d*x*/d*t*=0. Because the ODE form of Schlögl's model is a cubic, there can be one, two or three steady states for a given set of parameters. Note that the cubic term comes from the trimolecular portion of the reaction: bistability is not possible in a one-dimensional model involving only uni- and bimolecular reactions.

We will begin with the simplest case. Recall that, in general, in order to keep the concentrations of A and B at constants *a* and *b*, there will be a net flux of A going into the system and B coming out of the system, or vice versa. However, there is a particular ratio of *a* and *b* at which there will be no need for an external agent to add or remove A and B molecules from the system to maintain their concentrations on average (i.e. ignoring microscopic fluctuations). The values of *a* and *b* that make this possible are called equilibrium concentrations. These concentrations lead the system to a unique and completely describable fate: the equilibrium steady state.

*Chemical detailed balance* occurs when , so that the forward and backward fluxes are balanced in each and every reaction in the system. Note that this is a stronger condition than mathematical detailed balance, which simply assumes that the total forward rate is equal to the total backward rate. Using the definitions in equation (3.1), the chemical detailed balance (or equilibrium steady state) condition for Schlögl's model is(3.3)This condition is derived later in equation (3.10) by assuming that there is no chemical potential difference between A and B and no free energy lost when exchanging an A for a B or vice versa.

The value of *x* for which chemical detailed balance occurs is denoted as *x*_{ess}. This is the *equilibrium steady state*. When a system is closed to outside influences, or in contact with material baths with a constant chemical potential, this is the unique concentration value of *x* to which the system will tend. In the CME setting, the steady-state probability distribution of the equilibrium steady state is a Poisson (or multi-Poisson) distribution, as predicted by Gibbs' theory for grand canonical systems (Gardiner 1985; Heuett & Qian 2006). The number of molecules at which the distribution peaks corresponds (when divided by volume size) to the steady-state concentration in the ODE model. For our example, the steady-state concentration is(3.4)and the steady-state probability distribution is(3.5)

When the parameter set does not satisfy equation (3.3), the system will reach a NESS. Unlike equilibrium conditions, multiple NESS may exist for a given set of parameters. Since it is not necessary for , an NESS is characterized by the reactions forming a cycle: the reversible reaction in equation (2.1) makes an X from A, and the reaction in equation (2.2) transforms an X to B. Then, an external agent has to constantly make A from B to keep the concentrations of A and B constant. Hence, there is a chemical reaction cycle A→X→B→A. In the stochastic model in §3.3, we shall quantify this reaction cycle by a cycle flux.

We now categorize the types of behaviours that may occur under non-equilibrium conditions through bifurcation analysis, including a ‘bifurcation’ which can only be observed through the CME. Figure 1 shows plots of the ODE, d*x*/d*t* versus *x*, for an example of chemical equilibrium and for a bistable set of the parameters. When there is only one real root, the steady-state behaviour of the stochastic system will match that of the deterministic system (Mansour *et al.* 1981). However, as the parameters move away from the chemical detailed balance condition, the shape of the probability steady-state function will deform from its Poissonian shape.

In a biological setting, it is more likely that the concentrations of the pump parameters (A and B) could differ in varying situations, while the rate constants, which are inherent to the types of molecules involved, would stay the same. Thus, we assume that only the concentrations of A and B will change, and locate curves in the *ab*-plane at which the bifurcation to bistability occurs. The bifurcation point can be found through the discriminant of equation (2.3),(3.6)When *Δ*>0, there is only one steady state. When *Δ*<0, there are three steady states: two stable and one unstable. The curves for which *Δ*=0 in the *ab*-plane are shown in figure 2. Nicolis and Turner developed a method based on generating functions to analyse fluctuations in the steady state and have shown that, in the case of Schlögl's model, the variance diverges at the bifurcation point, *Δ*=0 (Nicolis & Turner 1977).

Under bistable conditions, the stochastic and deterministic models yield different predictions (an example of Keizer's paradox, see Vellela & Qian 2007). In the deterministic model, the system will tend towards one of the two stable fixed points, depending on the initial condition. The stochastic model also predicts that the system will quickly relax towards one stable point (FCA), but randomly switch between the two FCA on a time scale related to the system size.

Plotting the steady-state solution to the master equation (equation (2.8)), we observe that the two peaks can exchange their relative heights under changes in *a*, *b* or *V* individually (i.e. changing one parameter while keeping the other two constant). Figures 3 and 4 illustrate this switch under a change in the pump parameter *b*, and the volume size *V*, respectively. The volume size is a critical part of this bifurcation: the value of *a* or *b* at which the peaks have equal height depends on the volume size. Thus, this change cannot be observed through studying the ODE model.1^{,}2

It has been shown that under multistable conditions, it is impossible to derive a Fokker–Planck equation that correctly describes the behaviour predicted by the master equation (Horsthemke *et al.* 1977). The relative heights of the two peaks in the steady state are incorrectly approximated in the continuous limit by diffusion approximation. Since these heights dictate the relative stability of the FCA, and the height differences are amplified in the limit of large system size, accuracy of the probability steady state is critical.

None of the continuous models are able to accurately describe the switching behaviour for small volume systems, as we illustrate in §4.3. As of yet, there is no experimental data to confirm the predictions from any of these models. However, simulations through molecular dynamics (Baras & Mansour 1997) have shown agreement with the master equation. Thus, when studying non-equilibrium processes occurring on the microscopic and mesoscopic scales, the CME is the most valid mathematical model.

### 3.2 Deterministic chemical potential and entropy production

In this section, we attempt to characterize the steady-state behaviour using concepts from physical chemistry. In each of the reactions there is a *chemical potential difference*, between the Gibbs free energy stored in the system before the reaction occurs versus after it has occurred. The first reaction in the Schlögl model changes an A molecule into an X molecule. By definition (Beard & Qian 2008), the amount of chemical potential energy it takes to change one A into one X is(3.7)where *k*_{B} is Boltzmann's constant and *T* is the temperature (constant throughout our theory). Likewise, the Gibbs free energy used to change one B into one X is(3.8)The overall result of this reaction system is creating B molecules out of A molecules (or vice versa). Combining the two chemical potentials above, we see that the amount of potential energy lost in exchanging one A for one B is(3.9)where we have used *Δ**μ*_{XB}=−*Δ*μ_{BX}. We rename this expression *Δ**G* because it represents the change in the Gibbs free energy which occurs through one ‘forward’ cycle of the reaction (i.e. producing one B from one A).

The chemical potential difference depends on the values of the rate constants (enthalpic) and the amounts of substrate and product (entropic). For a special value of the parameters, namely(3.10)the chemical potential difference between A and B is zero. This is the same condition as in equation (3.3) and the system will tend towards the equilibrium steady state discussed in §3.1. In this state, there is no preference for producing B from A or vice versa because there is no energy difference between them.

Another important physical quantity is the *entropy production rate*. Entropy is created when an amount of useful energy is converted to heat in an isothermal process. In order to keep the concentrations of A and B constant, an external agent has to constantly remove the B, convert B back to *A,* and put A back in the system. All this work becomes heat in the steady-state reaction.

If chemical energy is constantly pumped into the system, a NESS will be reached in which the entropy is produced at a constant rate (related to the amount of work being done on the system) converted to heat and dissipated. The amount of entropy produced in turning one A into one B is equal to the free energy (differing by a constant factor of *T*, the temperature in Kelvin; equation (3.9)). Thus, the entropy production rate will be this constant multiplied by the rate of producing B from A. This rate is called the *steady-state flux*(3.11)where is the left, middle or right fixed point of equation (2.3) for *i*=−, 0 and +, respectively. The two differences in equation (3.11) represent the net amount of A being used, and the amount of B being produced, respectively, per unit time. These two rates will be equal in the steady state. Thus, we have the entropy production rate in the *i*th steady state as(3.12)Note that when the system is in the equilibrium steady state, there is only one entropy production rate and it is zero.

For all other cases, there may be multiple entropy production rates, and they will always be positive. To show this, we rewrite the entropy production rate in terms of each reaction, and use both definitions of the steady-state flux in equation (3.11):(3.13)(3.14)Note that a term such as ln(*x*/*y*)(*x*−*y*) is always greater than or equal to zero, since (*x*−*y*)>0 implies ln(*x*/*y*)>0 and (*x*−*y*)<0 implies . Thus, we know that *ξ*_{i}≥0 in all cases.

### 3.3 Stochastic entropy production and cycle flux

The entropy production rate is an important characteristic of the NESS. The entropy production rates can be different for each steady state in the deterministic model of a bistable system. In the stochastic model, there is a unique steady-state distribution, even in the case of bistability. Likewise, there is a unique entropy production rate, even when the system is in a NESS. In this section, we describe how the stochastic entropy production rate is defined and how it relates to the deterministic rate from above.

Let us ‘visualize’ the stochastic dynamics of the chemical reaction system, one reaction at a time. Every time a reaction occurs, there is an associated amount of chemical energy loss (to the solvent via heat) or gain (from the solvent via heat). Since the occurrences of the reactions are stochastic, the system's net heat output fluctuates. When the system reaches its stationary state, this heat output is *not* constant, but continues to increase on average (i.e. it increases at a constant rate). Thus, there is a mean stationary entropy production rate.

The entropy production rate for the stochastic model is calculated through a double sum: first over each possible state, then over each reaction (Luo *et al.* 2002). Thus, we must first split up the birth and death rates from equation (2.6) into the transition rates corresponding to each reaction. We define(3.15)(3.16)(3.17)and(3.18)so that and . We may then define the stochastic forward and backward steady-state fluxes analogous to those in equation (3.1):(3.19)where is the distribution given in equation (2.8).

The stochastic entropy production rate (as defined in Luo *et al.* (2002)) is(3.20)which represents summing the product of flux and chemical potential over all the states in the probability steady-state distribution for each reaction separately.

This definition can be rewritten in terms of the product of the Gibbs free energy change, *ΔG*, and the cycle flux, a quantity analogous to the steady-state flux from equation (3.11). Figure 5 illustrates how the rates from equation (3.15) affect the movement on the Markov chain. The net forward movement is equal to . This is defined as the *forward cycle flux*. Likewise, the net movement backward is equal to and is the *backward cycle flux*. These fluxes represent the rate at which the cycle turns over in the forward (i.e. A→B) and backward (i.e. B→A) directions in a given state *n*.

The forward and backward cycle fluxes must be equal when the system is in a steady state. Thus, the double sum in equation (3.20) can be rewritten as a single sum using either cycle flux, and the log terms can be combined,(3.21)Plugging in the definitions of the single reaction fluxes, the log term becomes the Gibbs free energy from equation (3.9),(3.22)This term represents the entropy production of one forward completion of cycle *n* (which is actually independent of *n*), since the Gibbs free energy is converted into entropy. The stochastic entropy production rate then becomes(3.23)

In this form, we see that the stochastic entropy production rate is intimately related to the entropy production rates defined in equation (3.12). Note that when we make the substitution *x*=*n*/*V*, the birth and death rates in equation (3.15) are related to the fluxes in equation (3.1) by(3.24)The stochastic entropy production rate can then be written in terms of the deterministic fluxes as(3.25)Assuming that the probability steady state is dominated by the values at its two peaks, we can approximate the entropy production rate by the sum of the two terms corresponding to those states:(3.26)where the *n*_{±} represent the integers which approximate the the closest and we have used the approximation around these states. Using the definition of the deterministic entropy production rate from equation (3.12), we have(3.27)which states that the stochastic mean entropy production rate is approximately the sum of the deterministic entropy production rates of the stable equilibria, weighted by their probabilities, times the system volume. Thus, when all other parameters are held constant, the stochastic entropy production rate will grow linearly with volume size.

As seen in the plots of *p*^{ss} in figures 3 and 4, even when there are two peaks in the probability steady state, one of the peaks is usually much more dominant (note the log scale). Owing to this, the stochastic entropy production rate will follow the deterministic entropy production rate of whichever equilibrium is more stable. Figure 6 shows a plot of the possible entropy productions rates (the *ξ*_{−,0,+} and *ξ*_{CME}) over a range of volumes. Only for a small range of the parameters will *ξ*_{CME} not match one of ; this is the range over which the peaks are comparable in height.

A topic of current interest is whether a general principle involving the entropy production rate can be stated (Callen 1957; Glansdorff *et al.* 1974; Keizer & Fox 1974; Jaynes 1986; Qian 2002; Dewar 2005; Martyusheva & Seleznev 2006): is the rate of entropy production always maximized, minimized or neither? Schlögl's model illustrates how the entropy production of a system in a NESS is not necessarily minimized (Callen 1957; Jaynes 1986) or maximized (Dewar 2005; Martyusheva & Seleznev 2006) with respect to the three choices of the deterministic entropy production rate. We see that the more stable of the two stable equilibria can have either a high or a low entropy production rate, as is discussed in Landauer (1975) and Andresen *et al.* (1984), and that the stochastic entropy production rate mimics this value. However, we note that the middle value of the entropy production rate, which corresponds to the unstable equilibrium, is not captured by the stochastic model. For this example, the stochastic entropy production rate will be either the maximum or minimum of the three choices, since the value of *ξ*_{0} is between *ξ*_{−} and *ξ*_{+}. It is not known if this is true for the systems in general.

Han and Wang have shown that the lower entropy production rates are indicative of less dissipation cost in a complex biochemical network, leading to a network that is more robust to perturbations. Thus, minimizing the entropy production may be a design principle in the biological systems (Han & Wang 2007, 2008) and/or an indication of stability (or robustness) of a given FCA. In our model, one measure of robustness is the difference between the peaks of the probability steady state and the trough between them. We refer to this difference as the barrier height, and define(3.28)where *n*_{0} corresponds to the position of the trough, *n*_{0}=*Vx*_{0}. By varying *a* and keeping the other parameters constant, we can observe how *ξ*_{CME} changes with respect to the quantities . Figure 7 illustrates this. The dark lines show the relationship between the stochastic entropy production and barrier height for the more stable FCA only. The negative slopes of these dark lines do suggest that the entropy production decreases as the more stable FCA becomes increasingly robust.

### 3.4 Bistable system as a two-state Markov process

When the system is bistable, there exists a non-zero possibility of jumping from one stable state to another in the stochastic model. When the system starts near one of the functional reaction system states (FCA), it will fluctuate around it for a long time. However, since stochastic fluctuations are always present, eventually they will push the system far enough away from one FCA so that it reaches the point corresponding to the unstable fixed point of the ODE model. Then the system will relax towards the other FCA and remain there for a long time.

This behaviour motivates a new method for viewing bistable systems. The infinite-state Markov chain can be compressed into a model with only two states: the two stable FCA. Movement near the FCA can be modelled with simple Gaussian processes (Keizer 1979), and jumping between the FCA can be predicted using the two-state Markov chain. This simplified view of the system will capture the behaviours on both time scales, and allow for a new, simpler method of simulation.

The key behaviour of a bistable system is the transition between the FCA. The transition rates from *x*_{−} to *x*_{+} have been approximated asymptotically by Hänggi.3 He showed that the rates depend exponentially on the volume size, so that as volume increases, the likelihood of a transition between the FCA decreases exponentially. These rates were also derived through a different method in (Hinch & Chapman 2005). The expression for the transition rate from *x*_{−} to *x*_{+} is (Hänggi *et al.* 1984; Hinch & Chapman 2005)(3.29)and the transition rate from *x*_{+} to *x*_{−} is(3.30)In these expressions, *Δu*_{±}=*u*(0)−*u*(*x*_{±}) where and is the second derivative of *u*(*x*) evaluated at *x*_{i}. The functions and from equation (3.1) are the total continuous forward and backward fluxes and the volume is represented as 1/*ϵ*.

We may think of the long-term behaviour as a two-state Markov chain whose infinitesimal matrix contains the transition rates(3.31)The vector * P* has two components:

*P*

_{−}and

*P*

_{+}are the probabilities of being near the FCA

*n*

_{−}and

*n*

_{+}, respectively. The process jumps between the two states, with a waiting time in each state whose average is inversely proportional to the transition rate into that state. The long-run probability of finding the process in the

*n*

_{+}state is

*r*

_{+}/(

*r*

_{+}+

*r*

_{−}), and the probability of finding it in the

*n*

_{−}state is

*r*

_{−}/(

*r*

_{+}+

*r*

_{−}).

The full CME model can also be represented as a matrix equation (see §4.1). To see how the system in equation (3.31) is related to the original Markov chain formulation, we examine the eigenvalues and eigenvectors of the matrix *A*. Because the columns of *A* sum to 0, the first eigenvalue is *λ*_{0}=0. The second eigenvalue is the sum of the transition rates, *λ*_{1}=−(*r*_{+}+*r*_{−}). These are also the first two eigenvalues of the matrix *Q* which represent the full chain (see §4.1). As discussed in the §4, these eigenvalues are the most important for the full system because they represent the long-term behaviour.

The right eigenvectors of the matrix *A* (normalized so that the sum of the terms equals one) are(3.32)and the left eigenvectors are(3.33)These four eigenvectors are connected to the first four eigenvectors of the matrix *Q* in equation (4.2), which represents the full system. Figure 8 shows an example of the eigenvectors **v**_{0}, **v**_{1}, **w**_{0} and **w**_{1} for the matrix *Q*.

The right eigenvectors of *Q*, **v**_{i}, represent contributions to the solution vector (see equation (4.3) in §4.1). The two right eigenvectors of *A*, **v**_{i} contain a shortened version of the main contributors, **v**_{0} and **v**_{1}. The right eigenvector **v**_{0} represents the probability steady-state distribution. The two values in the eigenvector **v**_{0} represent the proportion of probability around each FCA in the steady state. Thus, if we sum the first *n*_{0} entries of **v**_{0}, they will be equal to *r*_{+}/(*r*_{+}+*r*_{−}) (the first entry of **v**_{0}) and the sum of the rest of the entries will equal *r*_{−}/(*r*_{+}+*r*_{−}) (the second entry of **v**_{0}). Likewise, if we sum the first *n*_{0} entries of **v**_{1}, they will be equal to 1/2 and the sum of the rest will be −1/2, the two entries of **v**_{1}.

The vectors **w**_{0} and **w**_{0} are both entirely constant, since the equations * wQ*=0 and

*=0 are both satisfied by vectors of ones. The eigenvector*

**w**A

**w**_{1}comprises two constant parts and changes values quickly near the unstable point

*n*

_{0}. The ratio of these two constants is the same as the ratio of the two entries of

**w**_{1}. Thus, we see that the two-state process based on the transition rates between the FCA shares key characteristics of the original stochastic model.

Assuming that the system acts as a two-state process with the Gaussian noise around each FCA, we can model the probability distribution as the sum of two unnormalized Gaussian distributions. From Hinch & Chapman (2005), the zeroth-order asymptotic approximation to the steady-state probability distribution is(3.34)where the constant *A* can change to match the function properly near each FCA. The unnormalized Gaussian that best matches this function near *n*_{−} is4(3.35)and near *n*_{+}, it is(3.36)

We might consider the ratio of finding the two-state process in one state versus another. By integrating the approximations in equations (3.35) and (3.36) on the left- and right-hand sides of *n*_{0} (and using the identity ), we can approximate this ratio by(3.37)If we instead use the asymptotic equations for the transition rates themselves (equations (3.29) and (3.30)), we find(3.38)which is the same result. Thus, we conclude that, asymptotically, the steady-state behaviour is well approximated by the Gaussian noise around both FCA, with the process in equation (3.31) dictating the jumps between the two.

When *V* is large, a realistic system with any external noise will almost surely be found in only one of two FCA. Which state it ends up in depends solely on the parameter values (i.e. whether *u*(*x*_{+})>*u*(*x*_{−}) or vice versa) and not on the initial condition! Deterministically, which state the system ends up in depends solely on the initial condition. This is a significant contribution from the stochastic model and another example of the Keizer's paradox (Vellela & Qian 2007). Although there is a result from Kurtz (1976, 1978) which states that the solutions of the two types of models must match to a negligible probability in the thermodynamic limit for all finite time, there is no such statement when lim_{V→∞} and lim_{t→∞} are switched. Here, we observe again the disagreement when we assume the system has reached a steady state before taking the thermodynamic limit.

When working with a small volume system, jumping from one FCA to another is an event that may happen on a reasonable time scale. The system will reach its steady state (i.e. infinite time limit) without reaching the infinite volume limit. Thus, an understanding of the correct steady-state behaviour is crucial. In the §4, we derive a numerical method to examine small volume behaviours, and compare the asymptotic results with direct results from the master equation.

## 4. Numerical methods

In this section, we investigate a numerical method for analysing the Schlögl model. We show that Matlab is capable of accurately computing the eigenvalues and eigenvectors of a truncated linear system that represents the master equation. The exponentially small eigenvalue is well approximated for appropriately large truncations. This method not only produces the probability steady-state distribution; it also yields the full, time-dependent solution of the truncated system.

### 4.1 Calculating the solution

The stochastic model is represented by the infinite linear system(4.1)where , the vector of probabilities for each state. The matrix,(4.2)is a stochastic matrix, which has all real eigenvalues, the largest being *λ*_{0}=0.

Although all eigenvalues *λ*_{i} with *i*≥1 are real and negative, the second largest eigenvalue (or the eigenvalue with the second smallest magnitude), *λ*_{1}, differs from the others. This eigenvalue is exponentially small, i.e. it decays exponentially as the volume size increases. The other eigenvalues are relatively stationary with the volume change, and are much larger in magnitude for large volume. This behaviour is related to the exponential decay seen in the transition rates (see equations (3.29) and (3.30)). In fact, the exponentially small eigenvalue represents the sum of the two transition rates and therefore dictates the slow time scale of the system.

The solution in terms of the eigenvalues *λ*_{i} and eigenvectors **v**_{i} of *Q* will be(4.3)Because all of the *λ*_{i} are real and negative, all the terms in this solution will decay to zero as time increases except the first term, corresponding to *λ*_{0}=0. The eigenvector for the zero eigenvalue **v**_{0} therefore represents the steady-state probability distribution. The second eigenvalue *λ*_{1} is related to the time scale of the transitions between the FCA, and the contribution of its eigenvector **v**_{1} to the solution will decay on a much slower time scale than the eigenvectors **v**_{i} for *i*>1.

Although this system is infinite in size, the most probable states are those near the two FCA. The states in which *n* is much larger than *x*_{+}*V* are very unlikely, and the stationary probability distribution tends to zero as *n* tends to infinity. Therefore, it is reasonable to cut the system off at some finite value *N* and analyse the solution under these conditions (see Vellela & Qian (2007) for verification of this method). Thus, we consider the finite system formed by setting *λ*_{N}=0, so that the system cannot move past *N*. Along with setting *μ*_{N+1}=0, the corresponding transition matrix will remain a stochastic matrix with the expected properties.

We formed the matrix *Q* using the `spdiags` command in Matlab, with the following parameter values (which lie in the bistable region of parameter space):(4.4)The function is shown in figure 1 in §3.1. There are two stable steady states: *x*_{−}=0.0935 and *x*_{+}=3.7024i and one unstable steady state: *x*_{0}=1.2041.

The matrix was truncated at a value *N*, past the right fixed point *x*_{+}. The choice of *N* was based on a study of the convergence rate of the exponentially small eigenvalue for different choices of *N*. As in Vellela & Qian (2007), the eigenvalue quickly converges for truncation values past *N*=*x*_{0}*V*. The truncation value of *N*=2*x*_{+}*V* is sufficient, and was used for the rest of the calculations below. Eigenvalues and eigenvectors were calculated using the command `eig(full(Q))`.

### 4.2 Accuracy of the eigenvalues

The time-dependent solution can be expressed as in equation (4.3) as long as the calculated eigenvalues and eigenvectors are accurate. To test the accuracy of Matlab's results, we calculate the absolute residual vector,(4.5)and look at the *L*_{2} norm. Table 1 gives a list of the residuals for the first five eigenvalue/eigenvector pairs using a volume size of *V*=100 and the parameter values in equation (4.4). The absolute residual is of the order of 10^{−11} for all eigenvalues when *V*=100.

Other values of *V* were also examined, and it was found that the absolute residual grows slowly and approximately linearly for all eigenvalues as *V* increases. Although this residual is of the same order for all eigenvalues, *λ*_{1} grows exponentially with volume size, while the others remain relatively constant. A study of the relative residual, i.e. the absolute residual divided by the magnitude of the eigenvalue, may serve as a better measure of the accuracy of the Matlab's results. Figure 9 shows a plot of the relative residual, ‖*r*_{i}‖_{2}/|*λ*_{i}| for increasing values of *V*. Similar to the eigenvalues themselves, this residual grows exponentially with *V* for *λ*_{1} and is constant for the other eigenvalues. Since the magnitude of *λ*_{0} is so small, its relative residual is expectedly high. However, this residual remains constant with increasing *V*. For large volumes, such as *V*≥100, the absolute residuals are still small, but the relative residual for *λ*_{1} becomes greater than one and the exponential growth is lost. Thus, we cannot trust the Matlab's results for *V* in this region.

To check the accuracy of the Matlab's eigenvector calculations, we compare the eigenvector **v**_{0} with , which can be calculated exactly using equation (2.8). Figure 10 is a plot of the actual probability steady state versus the eigenvector **v**_{0}, for *V*=100. The eigenvector is a good approximation for the values above machine epsilon (10^{−16}). However, as *V* increases, the vector **v**_{0} begins to diverge from the true probability steady state for values near the FCA (as is seen around *n*_{+} in figure 10). A study of the difference between these two vectors, , shows that the error in **v**_{0} grows exponentially with volume size as well. As in figure 9, for volumes greater than 100, the error fluctuates wildly and even the exponential growth is no longer observable.

The first two right and left eigenvectors, **v**_{0}, **v**_{1}, **w**_{0} and **w**_{1}, respectively, are shown in figure 8. Another check of the Matlab's accuracy is to verify the orthogonality condition,(4.6)The left eigenvectors were calculated by calculating the right eigenvectors of the transpose matrix with `eig(full(Q).′)`, and taking the transpose of the eigenvector output. Both **v**_{0}.**w**_{1} and **v**_{1}.**w**_{0} were observed to grow exponentially with increasing *V*. For *V*=40, the dot products are(4.7)and(4.8)and for *V*=100, they are(4.9)and(4.10)Again, for volume sizes greater than 100, the error becomes too large to display the exponential growth with *V*.

These tests confirm that the Matlab's calculations are accurate for a range of small *V*, but as the volume and the matrix size grow, so do errors in the eigenvalues and the eigenvectors. The first two eigenvalues and their corresponding left and right eigenvalues are the most difficult to calculate, and the error grows exponentially as *V* increases. In the §4.3, we compare the results for *λ*_{1} from Matlab with the asymptotic predictions from Hänggi *et al.* (1984) and Hinch & Chapman (2005). As we will see, asymptotic approximations have the opposite problem: they are accurate for large *V* and break down for small *V*.

### 4.3 Comparison with simulations and asymptotics

Another way to check the accuracy of the exponentially small eigenvalue *λ*_{1} is to compare the Matlab's *λ*_{1} with a value obtained from the Monte Carlo simulations and with the sum of the asymptotic transition rates from equations (3.29) and (3.30). From Hänggi *et al.* (1984), the relationships between *λ*_{1}, the transition rates d*P*_{±}/d*t* and the waiting times *T*_{±} is(4.11)

The waiting times were calculated from the simulations by keeping track of which FCA the system was near (i.e. which ‘well’ the system was in), and the amount of time it spent while near that FCA. By dividing the total time spent in each well by the number of times that well was visited, we obtained the average waiting time near each stable FCA. The two waiting times were used in equation (4.11) as *T*_{+} and *T*_{−} to calculate the exponentially small eigenvalue.

The definition of one ‘visit’ to a well must be explicitly stated. The state *n*_{0}=round(*x*_{0}*V*) was first used as the dividing line: if the system was greater than *x*_{0}*V*, it was defined to be in the right-hand well, and if it was less than *n*_{0}, the left-hand well. However, it was observed during simulations that the system does not necessarily fall immediately into one well or another after passing over this point. Occasionally, the simulation would return to *n*_{0} several times before heading towards one of the FCA. This behaviour caused the calculation of average waiting time to be underestimated (and *λ*_{1} to be overestimated), because the number of visits was increased by one each time the state *n*_{0} was reached.

Thus, when calculating the average waiting times, we assumed that the simulation had visited a well only when it had reached a particular point, closer to the FCA than just *n*_{0}±1. A study was done of the change in *λ*_{1} with respect to the definition of this point. The eigenvalue displayed a fast, exponential-like convergence as the boundary moved away from *n*_{0}. For the rest of the simulations, we defined the point at which the system had officially entered a well to correspond to the inflection point in the potential function generated by integrating the cubic ODE.

In figure 11, we see agreement between the simulation predictions, the asymptotic approximations and the Matlab's calculations, and the exponential decay of *λ*_{1} as volume increases. The values obtained through the Monte Carlo simulations and the Matlab calculations are fairly similar. Although there is a small gap between these two lines, the slopes are very close, which indicates the same rate of predicted exponential decay with volume size.

However, the asymptotic predictions are much smaller than the other two, and the slope of the line is noticeably different. This is because the asymptotic equations are only valid for large volume systems. The simulations and the Matlab calculations were carried out for relatively small volumes, because of computing restrictions in both cases.

An attempt to unify the asymptotic results with those of the Matlab and the Monte Carlo was made. The slopes of the lines seen in figure 11 are related to the Δ*u*_{±} from equations (3.29) and (3.30). The constants Δ*u*_{±} are the asymptotic (continuous space) approximations to the barrier heights, from equation (3.28). Although it is too computationally intense to compute *λ*_{1} for large *V* using the `eig` command or a simulation, it is possible to compute from the detailed balance equation for larger *V* using Matlab. Thus, we can study the change in with volume. Figure 12 shows the convergence of the barrier heights to the asymptotic values Δ*u*_{±} as the system volume increases. We see the asymptotic prediction for the height differences in the probability steady state become accurate for volumes of the order of 100 and higher. Although the asymptotic line for *λ*_{1} in figure 11 is much lower than the other two lines, figure 12 tells us that eventually these lines should all converge as *V* increases further.

We can obtain the Matlab and the Monte Carlo data for values of *V* in the approximate range 5≤*V*≤40, and have shown that these two methods yield similar predictions for *λ*_{1}. We have also shown the range of *V* on which the asymptotic prediction becomes realistic. Unfortunately, these two ranges of volume size do not overlap. Computational limitations do not allow us to confirm the Matlab's predictions with those made in Hänggi *et al.* (1984) and Hinch & Chapman (2005) or vice versa. However, a study of this kind can suggest which method is the best when working with a model of a real biochemical system, in which the volume size is known.

## 5. Conclusion and future directions

Schlögl's model is a simple example of a chemical reaction system that exhibits bistability. Bistable behaviour can be found in many biological networks, including the heart models (Hinch & Chapman 2005), the visual perception (Moreno-Bote *et al.* 2007), the gene networks (Paulsson 2005), etc. Owing to the ubiquity of switching behaviour, it is important to have comprehensive mathematical models of the bistable chemical reaction systems. In this paper, we used Schlögl's model as an example to study the differences between the deterministic and the stochastic modelling techniques applied to bistability, and examine this behaviour through a thermodynamic perspective.

For a finite volume system, the steady-state behaviours predicted by the deterministic and the stochastic models can be drastically different. When the system is bistable, the behaviour in the deterministic model depends entirely on the initial condition, i.e. on which side of the unstable steady state the system started. On the contrary, the steady-state behaviour in the stochastic model is completely independent of the initial condition.

In the long term, the stochastic model predicts that the system spends almost all its time in the two FCA (states that correspond to the stable states of the ODE model), and the proportion of time spent in each is dictated by the ratio of the transition rates between them. Because this ratio increases exponentially with volume, as the volume size increases, the proportion of time the system spends in the more stable FCA increases as well. Thus, for a large but finite volume size, the stochastic models predicts that, independent of the starting point, the system will end up in the more stable FCA and spend most of its time there.

The parameters determine which of the FCA is more stable. Asymptotic analysis, which assumes a large volume limit, predicts that one of the FCA will be more stable for all values of *V*. Although this is true for large *V*, we have shown that for small *V*, the FCA may exchange relative stability with the changes solely in *V*. This observation highlights the importance of using a CME model when addressing the intracellular reaction networks. The stability exchange between the FCA has been observed with changes in the pump parameter *b* and the system volume. In the future, a full bifurcation diagram could be drawn so that the precise effects of *a*, *b* and *V* on the system's stability are known.

When the system is bistable, there are multiple entropy production rates that come from the ODE model, one for each steady state. After calculating this rate in the master equation model, we see that the stochastic entropy production rate is a combination of the rates in each possible state, weighted by the probability of being in each of those states. Thus, the stochastic rate ‘follows’ that of the deterministic system for the more stable FCA (which is the most likely state of the system). This shows that the entropy production rate is neither maximized nor minimized for all NESSs. We conclude, as in (Landauer 1975), that the entropy production rate alone is unable to fully describe stability in the steady state. However, this example does not disprove the conjecture that the entropy production rate will take on some critical value (i.e. a minimum or a maximum) in the NESS. A similar study of a more complicated model, such as a tristable model, may also debunk this idea.

The long-term behaviour of the bistable system is determined by the transition rates between the FCA. For large volumes, the system can be simplified by a two-state Markov chain in which all states to the left of the unstable fixed point are lumped into one state, and all states to the right are considered as another. The transition rates are the entries of the four-by-four stochastic matrix modelling the jumping between the FCA, and while near one FCA, we can model the short time-scale behaviour with Gaussian noise.

The sum of these two transition rates is equal to the second largest (or the second smallest in magnitude) eigenvalue in the stochastic matrix that represents the full Markov chain. This eigenvalue can be calculated with some accuracy using the `eig` command in Matlab. For small volume systems, the Matlab result is more accurate than the asymptotic prediction, and it is faster than calculating the eigenvalue through the Monte Carlo simulations. However, a reliable method has not yet been found to compute this eigenvalue for all ranges of volume. Such a method will be crucial in predicting the time scales of a bistable chemical reaction system for mesoscopic volume sizes.

## Acknowledgments

We thank Profs Ping Ao, Elliot Elson and Annie Lemarchand for their helpful discussions and Joshua Goldwyn and Pei-Zhe Shi for reading the manuscripts. The work is partially supported by the NIH grant no. GM068610.

## Footnotes

↵For example, this bifurcation does not correspond to the point at which the potential function has two minima of equal depth.

↵In a more strict sense, the relative stability of the two peaks depends not only on the heights of the peak, but also on the curvature. In thermodynamic terms, the former is enthopic and the latter is entropic contribution to the free energy. For large system size, the dominant effect is from the height difference, which is exponentially related to the volume

*V*. The curvature converges when*V*→∞.↵Although the transition rates in Hanggi's model match the rates predicted by the master equation in the limit of large system size, it is important to note that the solution to Hanggi's model does not match the solution to the CME at any point other than in the steady state, and that his rates are not a good approximation in small volume systems (see §4.3).

↵This approximation was found by expanding the logarithm of the function in equation (3.34) about

*n*_{±}and matching the first three terms to fit a Gaussian, noting that*J*^{+}=*J*^{−}at the fixed points.↵For example, using the parameter values

*k*1=3,*k*2=0.6,*k*3=0.25,*k*4=2.95,*a*=1, and*b*=1.305, the two integrals have different signs.- Received November 5, 2008.
- Accepted November 25, 2008.

- © 2008 The Royal Society