Hello,

I am studying a basic model that will aid in understanding more complex models. The basic model is the Michaelis-Menten Reaction with Substrate Input with No Product. The enzymatic reaction scheme is below:

Null -> S + E <-> SE -> E,

with production rate ‘Q’, complex association rate ‘a’, complex dissociation rate ‘d’, and catalytic rate ‘k’.

I wrote down the Reaction Rate equations (ODEs) for the system and solved for the steady-state concentration of S. It turns out that it is defined only when Q < V_{max,E}. There is a “protein catastrophe” for large enough Q.

This is an open system that does not satisfy the property of detailed balance. I tried investigating with an analytical treatment of the Chemical Master Equation. I also tried Moment Closure Approximation methods. After much effort and literature survey, I found it is highly unlikely that one can derive the stationary distribution for S, let alone its limiting distribution.

So I sampled trajectories of S from the (exact) CME using the Gillespie algorithm. I began the stochastic simulation at the steady-state concentrations given by numerical integration of the ODEs. I chose parameter values such that the system was very close to the protein catastrophe (i.e. large enough Q to yield S ~ 10^7). However, even after running the simulation for a very long time, I found that both mean(S) and variance(S) had not settled to a steady-state. Instead, they were fluctuating with respect to simulation time.

My questions:

1. Why could this be happening for the Mean? I thought that the large species numbers in the system would remedy any discrepancy between the deterministic average and the stochastic mean.
2. Why for the Variance?

You see, this poses some issues because my end goal is to calculate the Coefficient of Variation for S, as a function of Q.

Thanks for any help.