.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/kinetics/networks/plot_06_gillespie_stochastic_simulation.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_kinetics_networks_plot_06_gillespie_stochastic_simulation.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_kinetics_networks_plot_06_gillespie_stochastic_simulation.py:


Gillespie's stochastic simulation of reaction kinetics
=========================================================

Gillespie (1976-1977) showed how to sample *exact* trajectories of a
reaction network treated as random, discrete molecular events: draw an
exponential waiting time with rate equal to the total propensity, then
pick which reaction fires in proportion to its propensity.
:func:`~chemistrykit.kinetics.systems.stochastic.gillespie_ssa` implements
this direct method. With thousands of molecules the trajectories hug the
deterministic rate-equation solution; with a handful they are visibly
noisy, and the average over many runs recovers the deterministic curve.

.. GENERATED FROM PYTHON SOURCE LINES 16-49



.. image-sg:: /api/gallery/kinetics/networks/images/sphx_glr_plot_06_gillespie_stochastic_simulation_001.png
   :alt: A -> B -> C with N = 20 molecules (5 sample paths), A -> B -> C with N = 2000 molecules (5 sample paths)
   :srcset: /api/gallery/kinetics/networks/images/sphx_glr_plot_06_gillespie_stochastic_simulation_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    N=20: max |mean SSA - deterministic| = 0.0319
    N=2000: max |mean SSA - deterministic| = 0.0013






|

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from chemistrykit.kinetics.systems.networks import consecutive_analytic
    from chemistrykit.kinetics.systems.stochastic import gillespie_ssa

    # A -> B (c1), B -> C (c2): counts of molecules, not concentrations.
    stoich = [[-1, 0], [1, -1], [0, 1]]
    orders = [[1, 0], [0, 1], [0, 0]]
    c1, c2 = 1.0, 0.3
    species = ("A", "B", "C")
    t_grid = np.linspace(0.0, 15.0, 301)
    rng = np.random.default_rng(1976)

    fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
    for ax, N in zip(axes, (20, 2000)):
        _, B_det, _ = consecutive_analytic(1.0, c1, c2, t_grid)
        runs = []
        for i in range(100):
            traj = gillespie_ssa(stoich, [c1, c2], orders, [N, 0, 0], t_max=15.0, species=species, seed=rng)
            runs.append(traj.sample(t_grid)[:, 1] / N)
            if i < 5:
                ax.step(traj.t, traj.count("B") / N, where="post", color="steelblue", alpha=0.5, linewidth=0.8)
        ax.plot(t_grid, np.mean(runs, axis=0), color="crimson", label="mean of 100 SSA runs")
        ax.plot(t_grid, B_det, "k--", label="deterministic rate equations")
        ax.set_xlabel("t")
        ax.set_ylabel("[B] / N_A(0)")
        ax.set_title(f"A -> B -> C with N = {N} molecules (5 sample paths)")
        ax.legend()
        print(f"N={N}: max |mean SSA - deterministic| = {np.max(np.abs(np.mean(runs, axis=0) - B_det)):.4f}")

    fig.tight_layout()
    plt.show()


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.671 seconds)


.. _sphx_glr_download_api_gallery_kinetics_networks_plot_06_gillespie_stochastic_simulation.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_06_gillespie_stochastic_simulation.ipynb <plot_06_gillespie_stochastic_simulation.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_06_gillespie_stochastic_simulation.py <plot_06_gillespie_stochastic_simulation.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_06_gillespie_stochastic_simulation.zip <plot_06_gillespie_stochastic_simulation.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_