The stochastic reactions for R*, described in the panels in Figure 1, were simulated using Gillespie’s method [19 ]. We will begin by describing the simplest cases (Figure 1A–C) where just two fates are possible for each state: either phosphorylation or arrestin binding. First, the mean transition rate constants μ(n) and ν(n) were assigned. Then, for each states, two pseudo-random numbers, r1(n) and r2(n), uniformly distributed between zero and unity, were used. The lifetime that the molecule remained in state R*·P(n) was simulated as −ln(r1(n))/{ν(n)+μ(n)}, and the state to which the molecule transitioned was simulated as phosphorylation if r2(n) < ν(n)/{ν(n)+μ(n)} and as arrestin binding otherwise. If a finite flash width was required, the time of photoisomerization was simulated as uniformly distributed over that flash width. These simulations were very fast, and 106 iterations could be performed in ~1 s on a laptop PC. For the three-state R* activity model, simulation of the full version (Figure 1D) would have required more complicated code, so we chose instead to simulate only the truncated linear configuration (Figure 1E), and in this case 106 iterations again took only ~1 s.