Comparison of the kinetics of different Markov models for ligand binding under varying conditions

We recently derived a Markov model for macromolecular ligand binding dynamics from few physical assumptions and showed that its stationary distribution is the grand canonical ensemble [J. W. R. Martini, M. Habeck, and M. Schlather, J. Math. Chem. 52, 665 (2014)]. The transition probabilities of the proposed Markov process define a particular Glauber dynamics and have some similarity to the Metropolis-Hastings algorithm. Here, we illustrate that this model is the stochastic analog of (pseudo) rate equations and the corresponding system of differential equations. Moreover, it can be viewed as a limiting case of general stochastic simulations of chemical kinetics. Thus, the model links stochastic and deterministic approaches as well as kinetics and equilibrium described by the grand canonical ensemble. We demonstrate that the family of transition matrices of our model, parameterized by temperature and ligand activity, generates ligand binding kinetics that respond to changes in these parameters in a qualitatively similar way as experimentally observed kinetics. In contrast, neither the Metropolis-Hastings algorithm nor the Glauber heat bath reflects changes in the external conditions correctly. Both converge rapidly to the stationary distribution, which is advantageous when the major interest is in the equilibrium state, but fail to describe the kinetics of ligand binding realistically. To simulate cellular processes that involve the reversible stochastic binding of multiple factors, our pseudo rate equation model should therefore be preferred to the Metropolis-Hastings algorithm and the Glauber heat bath, if the stationary distribution is not of only interest.


I. INTRODUCTION
The function of many biological macromolecules is modulated by interactions with other molecules.Examples include enzyme cofactors as well as ligands causing heterotropic allosteric effects such as 2,3-bisphosphoglycerate, which changes the affinity of oxygen to the binding sites of hemoglobin.Experimental studies of specific molecule-ligand systems and theoretical investigations into the mechanisms of ligand binding both help us to understand the nature of (bio)chemical processes.
][3] However, in addition to the equilibrium properties of the molecule-ligand system, one is also interested in questions related to the kinetics and statistics of the binding process: How fast is the equilibrium state reached?What is the effect of a change in an external parameter such as temperature or ligand concentration on the kinetics?How does the variance in the number of bound ligands change over time?a) Electronic mail: jmartin2@gwdg.deb) Electronic mail: mhabeck@gwdg.de   One way to study the kinetics of ligand binding is to use molecular dynamics (MD) simulations.MD provides highly detailed information but is typically computationally very demanding and mostly samples events that are irrelevant to the binding process.If the focus is on the binding events themselves, stochastic methods for chemical kinetics are to be preferred over deterministic descriptions (for a review on stochastic simulation of chemical kinetics see Gillespie  (2007)  4 ).The stochasticity of the system will be important if the number of reacting molecules is small and the system cannot be characterized adequately by an average behavior. 5n stochastic simulations of chemical kinetics, usually the entire system of reactants is considered and characterized by a state vector with x i, t denoting the number of molecules of species i at time t.The stochastic dynamics of the system is described by the chemical master equation defining time-dependent probabilities that the system adopts a particular state. 6,7Many algorithms such as the stochastic simulation algorithm (SSA 8 ) have been developed to simulate trajectories of the whole system.
In the following, we assume that the reactants are wellstirred and that the number of ligand molecules is much greater than the number of target molecules.This assumption does not necessarily justify the use of deterministic models, because the number of target molecules can still be small and thus the stochasticity of the system may be important.However, we can make the following simplification: Upon ligand binding, the number of available ligands will be reduced and the probability of another binding event will be slightly smaller.This means that the states of the target molecules are stochastically dependent, since the state of a molecule influences the behavior of the other molecules.This coupling is fully taken into account if the whole system is described by state vectors of shape Eq. ( 1) and the chemical master equation.However, if we assume that the number of ligands is much greater than the number of target molecules, the effect of a reduction in the number of ligands on the transition probabilities can be neglected.The target molecules are then stochastically independent to good approximation, which implies that we can represent the entire system by a single molecule.In addition to an enormous reduction of the sample space, this simplification has the advantage that all reactions are of first order (pseudo unimolecular), which facilitates the calculation of quantities such as expected values (see Gillespie (2007)  4 ).
Here, we compare the stochastic pseudo rate equation model 9 to other Markov chains that share the same stationary distribution and that are commonly used in Monte Carlo simulations.We show that the stochastic pseudo rate equation model bridges between the general stochastic simulation of the entire system, pseudo first order rate (or differential) equations and the stochastic description of equilibrium by the grand canonical ensemble of statistical mechanics.

II. EQUILIBRIUM STATE DESCRIBED BY THE GRAND CANONICAL ENSEMBLE
In the grand canonical ensemble, the probability of finding a target molecule in a particular microstate k = (k 1 , . . ., k n ) ∈ {0, 1} n is the ratio of the Gibbs factor and the partition function (or binding polynomial) 1,2 with Boltzmann factor and binding polynomial where G(k) is the free energy of the molecule in state k compared to a fixed, arbitrarily chosen reference state, R the Boltzmann constant, T the absolute temperature.Morever, λ denotes the thermodynamic activity of the ligand and |k| the number of occupied sites.The derivation of the grand canonical ensemble (Eq.( 2)) assumes that the binding or release of a ligand molecule does not change the temperature or the chemical activity of the environment. 2This corresponds to the assumption that the target molecules are stochastically independent.Therefore, whenever the grand canonical ensemble describes the equilibrium state of an ensemble of target molecules adequately, the time evolution of the whole reaction system should be described well enough by a Markov process modeling the dynamics of a single target molecule.Questions such as "how rapidly does the system respond to a change in the external conditions such as ligand concentration?," can be answered if this Markov process is known.Following this reasoning, we compare the properties of a Markov process that we derived from few physical assumptions 9 to other Glauber dynamics with similar structure.Moreover, we discuss the analogy to (pseudo) rate equations and compare the effects of changes of the kinetics under altered environmental conditions to experimental data on myoglobin. 10The comparison provides evidence that the pseudo rate equation model is a qualitatively correct description of the dynamics of ligand binding.Finally, we use the model to study the effect of cooperativity on the kinetics of ligand binding.

III. DISCRETE-TIME MODELS OF LIGAND BINDING BASED ON GLAUBER DYNAMICS
Let us first introduce the Markov process from which we derived the grand canonical ensemble and compare it to two other well-known Glauber dynamics. 11The term "Glauber dynamics" is used in different fields of probability theory and statistical physics for not completely coinciding dynamics.Here, we use this term for any Markov chain on {0, 1} n whose transition probability is the product of a proposal and a conditional acceptance probability where the proposal distribution is uniform over all neighboring states.Following Martinelli, 12 we call a specific dynamics defined below the "Glauber heat bath."The three models that will be compared share the same proposal probability which is nonzero only for transitions k → l involving a single binding event |k − l| = 1.The only difference is in their acceptance probability.

A. Metropolis-Hastings algorithm
The Metropolis-Hastings algorithm 13,14 uses the acceptance probability and is commonly used in statistical simulations to "draw" samples from its stationary distribution π(k) [Eq.( 2)].

B. Pseudo rate equation analog
Our Markov model for ligand binding dynamics resembles the Metropolis-Hastings algorithm. 9The major difference is that we assume that a transition from a current state k to a neighboring state l involves two stochastically independent events: The environment has to provide or take up a ligand molecule, and the target molecule has to bind or release that ligand molecule.Therefore, the total acceptance probability is given by in case a ligand molecule is taken up by the environment, and when a ligand molecule is released into the environment.The factor min{1, g(l)/g(k)} represents the capability of the molecule to undergo the transition, whereas θ 1 and θ 2 are the probabilities that the environment provides or takes up a ligand molecule, respectively; their ratio can be interpreted as ligand activity λ = θ 1 /θ 2 . 9Thus, the only difference to the Metropolis-Hastings algorithm is that the factor λ is not included in the min-operator, meaning that we decompose ligand binding or release into two independent events.

C. Glauber heat bath
Another well-known dynamics is the Glauber heat bath. 11,12,15,16After a candidate binding site that might change its state has been chosen randomly, the acceptance probability is given by

A. Discrete-time models
For a molecule with only a single binding site, g 1 = exp {−G(1)/RT } ≥ 1 and θ 2 = 1, the transition matrices of the three Markov models are By construction, the stationary distribution π of all three Markov chains is (π 0 , π 1 ) = ( 1 1+g 1 λ , g 1 λ 1+g 1 λ ); it is the left eigenvector of P i with eigenvalue 1 (πP i = π).The second eigenvalue ν i determines the speed of convergence of the Markov chains.
The second eigenvalue of P 3 is zero (ν 3 = 0), which implies that the Glauber heat bath converges within a single transition to its stationary distribution, independent of the ligand activity λ and starting distribution.This behavior does not seem realistic for chemical reaction systems.However, by studying the continuous-time version of the Markov chains it might be possible to "zoom" into the details of a single transition and observe more realistic dynamics.

B. Continuous-time models
Continuous-time Markov processes are commonly used when the chemical kinetics is described by a master equation.
Moreover, it is conceivable that the continuous-time models show less abrupt changes in the system's state.The transition probability matrix of the Metropolis-Hastings algorithm, for example, has a negative second eigenvalue leading to strong fluctuations as the system approaches the stationary distribution.These fluctuations are reduced in the continuoustime analog, because the transition time is random and not all molecules change their state simultaneously.We will also see that the continuous-time version of the Glauber heat bath no longer converges within a single time increment.
The transition rate matrices of the continuous-time processes are given by with I denoting the identity matrix.According to the theory of continuous-time Markov processes on discrete spaces, [17][18][19] the transition probabilities within a time interval t are given by the matrix exponential where P i (t) solves the master equation Ṗi = P i Q i .The continuous-time Markov processes (11) share with the corresponding original discrete-time processes (9) the stationary distribution, the expected time a molecule stays in a particular state, and the distribution across the other states.The general structure of the time-dependent 2 × 2 matrix of transition probabilities is 20 where we introduced the parameters α i to allow for an appropriate scaling of t.The reaction rates r i = 1 − ν i are determined by the second eigenvalue ν i of the Markov models and given by for the Metropolis-Hastings algorithm (assuming g 1 λ > 1), the pseudo rate equation model (assuming g 1 > 1), and the Glauber heat bath, respectively.If the molecule is empty initially, the probability of being occupied after time t is To see whether the models are physically realistic, we will compare the kinetics to experimental data and investigate whether the behavior under altered environmental conditions (changes in λ and T) is reproduced by the models (see Sec. IV D).

C. Rate equations
Let us now study the binding reaction M + L M L using rate equations where M and M L denote the free and occupied target molecule and L the ligand.We assume that the activity or concentration of the ligand does not change significantly, because the number of target molecules is much smaller than the number of ligands.Upon ligand binding, the change in the thermodynamic activity (or concentration) µ of the empty target molecule is μ = −κλµ; κ scales the reaction time (and potentially also incorporates an activation energy barrier).The dissociation reaction is μ = κ(µ 0 − µ)/g 1 with µ 0 denoting the activity of all target molecules .The net change in free molecule concentration is given by The time evolution given by P 2 (t) solves Eq. ( 14): If the probability to be in the free state is p 0 initially, this probability will change to under P 2 after an evolution time t.In case there is a large number of target molecules, p 0 (t) can be interpreted as the fraction of free molecules identified with µ(t)/µ 0 .If we let κ = α 2 , p 0 (t) indeed solves Eq. ( 14) for any initial value p 0 .The same holds for the fraction of occupied molecules [µ 0 − µ(t)]/µ 0 .In contrast, the time evolution given by P 1 (t) and P 3 (t) is incompatible with Eq. ( 14).][23] Note that activation energies have not been considered so far.This can be remedied by substituting the time scales α i with temperature-dependent scales α i (T), because the reaction time depends on the probability to overcome the activation energy barrier, which itself will be a function of the temperature.

D. Application to ligand binding of myoglobin
We will now compare the continuous-time models to the ligand binding kinetics of myoglobin.8][29] Here, we only focus on the binding reaction and its dependence on external parameters and compare the qualitative behavior of the presented models to process IV according to Austin et al., 10 which exhibits singleexponential kinetics.In their figures 4 to 9, Austin et al. show the number of free molecules N(t) as a function of time, which corresponds to the first entry of the first row of P i (t).We assume g 1 λ ≫ 1 thereby approximating the stationary distribution by (π 0 , π 1 ) = (0, 1), because nearly all myoglobin molecules seem to be occupied in the equilibrium state.It is the decay rate r i [Eq.(13)] that determines the shape of the kinetics; the approximation (π 0 , π 1 ) = (0, 1) has only a marginal effect, e.g., in a logarithmic plot it determines the baseline of the reaction curve.The kinetics of process IV is then given by according to the Metropolis-Hastings, pseudo-rate and heat bath model, respectively.In the following, we will study the effect of changing the ligand activity or temperature.It is evident that the heat bath dynamics will not describe these changes adequately, because N 3 (t) is independent of ligand activity and binding energy.We will therefore mainly compare the Metropolis-Hastings and our pseudo rate equation process.

Changes in ligand activity
Let us first compare the effect of a change in ligand activity λ with the experimental curves.Inspection of Figure 6 in Austin et al. 10 highlights that changing λ by a factor of 100 enhances the binding process significantly.The Metropolis-Hastings process cannot explain this effect, because an increase in λ reduces the binding velocity according to P 1 (t).Therefore, only the pseudo rate equation model P 2 (t) seems suitable to describe the binding behavior of myoglobin.
Although Austin et al. 10 present many experimental data, it is not straightforward to estimate the parameters of our model, which would allow for a quantitative comparison.The reason is that multiple conditions (the type of solvent and ligand as well as concentrations and temperature) are varied simultaneously.However, the ligand activity and binding energy will both depend on the solvent and thereby have an influence on the binding kinetics. 30Other problems arise from the fact that the raw data are not available.We can therefore only show a qualitative comparison of our model with the findings of Austin et al.
We estimate the binding curves of process IV (Figure 6 of Austin et al. 10 ) at CO concentrations of 3 µM and 300 µM to follow According to our model, this effect implies that g 1 = 29/70λ where λ denotes the thermodynamic activity of CO in polyvinyl alcohol (PVA) at 3 µM concentration (see Fig. 1 for an illustration of the effect).

Changes in temperature
Austin et al. also study the effect of a change in temperature on myoglobin in glycerol-water (see their Figure 7).We assume that temperature changes mainly affect g 1 .The Gibbs free energy of the binding reaction is G = H − T S and favors ligand binding (G < 0).According to Eq. ( 3), g 1 = exp{−H/RT + S/R} ≥ 1.We further assume that ligand This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.76.223.157On: Mon, 13 Apr 2015 12:09:18 FIG. 2. Fraction of unoccupied molecules N (t) as given in Eq. ( 17) with q = 0.5, λ = 10 −5 , g 1 = 10 7 , and α 2 = 10 4 .
binding reduces the entropy S of the system which then implies H < 0.Then, g 1 will increase if T decreases.As a consequence, the binding rate λ + 1/g 1 will decrease according to our model.This behavior corresponds qualitatively to the curves shown in Figure 7 of Austin et al. 10

Low temperature and solid environment
The binding of myoglobin is comprised of several intermediate steps that cannot be described by a single exponential decay. 10Let us briefly outline how similar effects can be generated based on our process.Not all molecules may obey the time evolution 2 after photo-dissociation, because the ligand cannot diffuse away from the molecule due to the low temperature or solid environment.In this case, the ligand is immediately available for the corresponding fraction of target molecules q(T), where q(T) will increase monotonically in temperature T. These molecules will follow the law P 2 (t) |λ=1 such that overall The resulting curve is shown in Fig. 2; consider also the ligand binding of neuroglobin 31 for analogous curves.In the case of myoglobin which undergoes multiple intermediate transitions, the measured curve should be the sum of multiple exponential decays.

V. TWO BINDING SITES AND COOPERATIVITY
The state space of a molecule with two binding sites is {0, 1} 2 .The stationary distribution (2) can be parameterized by where g i are the binding constants (i.e., −RT ln g i are the free energy differences between states (0, 0), (1, 0), and (0, 0), (0, 1)) and w is a coupling term (i.e., −RT ln w is the free energy difference between (1, 1) and the sum of (0, 1) and (1, 0)).We will now study some implications of the pseudo-rate model for the role of cooperativity in the kinetics of ligand binding.To do so, we analyze three molecules with the same stationary distribution over macrostates |k| ∈ {0, 1, 2} but different interaction energies.Our reference molecule is , w M 1 ) = (10 6 , 10 6 , 10 −3 ) with binding polynomial 10 9 λ 2 + 2 • 10 6 λ + 1. Molecule M 1 shows negative cooperativity since w < 1. Molecules with the same equilibrium properties but different couplings w are obtained by solving g 1 + g 2 = g To study the effect of positive cooperativity, we choose an interaction constant of w = 10 3 > 1 resulting in the molecule 32 M 3 ≈ (1 999 999, 0.500 000 1, 10 3 ).
We are interested in how cooperativity affects the binding dynamics.To answer this question, we simulate Markov model (6) and, as an alternative approach, solve the deterministic rate equation ṗ = p Q for all three molecules where Q is the 4 × 4 rate matrix (recall that the deterministic approach equals the expectation of the stochastic model).As an observable, we choose the time evolution of the number of ligands bound to the target molecule, which in the deterministic approach is given by n We start with an empty molecule, i.e., in the deterministic approach, the initial distribution is p(0) = (1, 0, 0, 0).For comparison, we also simulate an ensemble of n independent target molecules stochastically using Markov model (6).Typical results are shown in Fig. 3 for systems comprised of 100 and 1000 target molecules.The observation that negative cooperativity leads to an acceleration of the binding kinetics was also confirmed by other simulated examples.Figure 3 also illustrates the correspondence between stochastic simulations and the deterministic rate equations.The average number of bound ligands fluctuates about the deterministic curve, showing that both agree better and better as the number of target molecules increases.

VI. SUMMARY AND OUTLOOK
In this article, we argued that under certain conditions the stochastic model of a reaction system can be reduced to a model for a single molecule, which is valid, for example, if the grand canonical ensemble describes the equilibrium state of the system of target molecules.This simplification reduces the state space of the Markov model drastically and saves computation times.The approach also allows for analytical calculations, because all reactions are pseudo unimolecular.Following this reasoning, we compared a Markov process that we found to be analogous to pseudo-rate equations to two other Glauber dynamics with similarly structured transition probability matrices.Using examples of ligands binding to myoglobin, we demonstrated that the kinetics of the Metropolis-Hastings algorithm and the Glauber heat bath are at odds with experimentally measured kinetics under varying external conditions, whereas the pseudo rate equation model shows the same qualitative behavior as the experiment.Additionally, we highlighted that the pseudo-rate equation process establishes a link between general stochastic simulations of whole reaction systems, rate equations, and the grand canonical ensemble.In particular, the model relates kinetic to equilibrium properties.For instance, N 2 (t) in Eq. (16) shows that altering the concentration λ 0 by a factor a has a larger effect, if the reference activity fulfills λ 0 ≫ 1 g 1 than if both quantities are approximately of equal size.Studies of systems with two binding sites illustrate that negative binding cooperativity may reduce the time until the equilibrium is reached and that in principle quantities such as interaction energies could be estimated from the kinetics of ligand binding and titration curves.In summary, our results indicate that the pseudo-rate equation model is more appropriate than the standard Metropolis-Hastings algorithm or the Glauber heat bath, if the kinetics itself is of interest and that it offers a simple alternative to a stochastic simulation of the entire system if the concentrations of the ligands are much higher than the concentrations of the macromolecules.
FIG.1.Effect of changing the ligand activity on N (t).The CO concentration ranges from 3 µM (black curve) to 300 µM (red curve).The dashed lines show N (t) at intermediate concentrations.