Algorithms

cayenne currently has 3 algorithms:

  1. Gillespie’s direct method (direct) (accurate, may be slow)
  2. Tau leaping method (tau_leaping) (approximate, faster, needs to be tuned)
  3. Adaptive tau leaping method (experimental, tau_adaptive) (approximate, faster, largely self-tuning)

Methods are described more in depth below.

Gillespie’s direct method (direct)

Implementation of Gillespie’s Direct method. This is an exact simulation algorithm that simulates each reaction step. This makes it slower than other methods, but it’s a good place to start.

class cayenne.algorithms.direct

Runs the Direct Stochastic Simulation Algorithm [1]

Parameters:
  • react_stoic ((ns, nr) ndarray) – A 2D array of the stoichiometric coefficients of the reactants. Reactions are columns and species are rows.
  • prod_stoic ((ns, nr) ndarray) – A 2D array of the stoichiometric coefficients of the products. Reactions are columns and species are rows.
  • init_state ((ns,) ndarray) – A 1D array representing the initial state of the system.
  • k_det ((nr,) ndarray) – A 1D array representing the deterministic rate constants of the system.
  • max_t (float) – The maximum simulation time to run the simulation for.
  • max_iter (int) – The maximum number of iterations to run the simulation for.
  • volume (float) – The volume of the reactor vessel which is important for second and higher order reactions. Defaults to 1 arbitrary units.
  • seed (int) – The seed for the numpy random generator used for the current run of the algorithm.
  • chem_flag (bool) – If True, divide by Na (Avogadro’s constant) while calculating stochastic rate constants. Defaults to False.
Returns:

  • t (ndarray) – Numpy array of the time points.

  • x (ndarray) – Numpy array of the states of the system at times in in t.

  • status (int) – Indicates the status of the simulation at exit.

    1 - Succesful completion, terminated when max_iter iterations reached.

    2 - Succesful completion, terminated when max_t crossed.

    3 - Succesful completion, terminated when all species went extinct.

    -1 - Failure, order greater than 3 detected.

    -2 - Failure, propensity zero without extinction.

References

[1]Gillespie, D.T., 1976. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22, 403–434. doi:10.1016/0021-9991(76)90041-3.

Tau leaping method (tau_leaping)

Implementation of the tau leaping algorithm. This is an approximate method that needs to be tuned to the system at hand (by modifying the time step given by the tau parameter). A default tau=0.1 is assumed by cayenne. This algorithm is approximate and faster than the Direct algorithm, but it must be used with caution. Smaller time steps make the simulation more accurate, but increase the code run time. Larger time steps make the simulations less accurate but speeds up code run time.

class cayenne.algorithms.tau_leaping

Runs the Tau Leaping Simulation Algorithm. Exits if negative population encountered.

Parameters:
  • react_stoic ((ns, nr) ndarray) – A 2D array of the stoichiometric coefficients of the reactants. Reactions are columns and species are rows.
  • prod_stoic ((ns, nr) ndarray) – A 2D array of the stoichiometric coefficients of the products. Reactions are columns and species are rows.
  • init_state ((ns,) ndarray) – A 1D array representing the initial state of the system.
  • k_det ((nr,) ndarray) – A 1D array representing the deterministic rate constants of the system.
  • tau (float) – The constant time step used to tau leaping.
  • max_t (float) – The maximum simulation time to run the simulation for.
  • volume (float) – The volume of the reactor vessel which is important for second and higher order reactions. Defaults to 1 arbitrary units.
  • seed (int) – The seed for the numpy random generator used for the current run of the algorithm.
  • chem_flag (bool) – If True, divide by Na (Avogadro’s constant) while calculating stochastic rate constants. Defaults to False.
Returns:

  • t (ndarray) – Numpy array of the time points.

  • x (ndarray) – Numpy array of the states of the system at times in in t.

  • status (int) – Indicates the status of the simulation at exit.

    1: Succesful completion, terminated when max_iter iterations reached.

    2: Succesful completion, terminated when max_t crossed.

    3: Succesful completion, terminated when all species went extinct.

    -1: Failure, order greater than 3 detected.

    -2: Failure, propensity zero without extinction.

    -3: Negative species count encountered.

Adaptive tau leaping method (experimental, tau_adaptive)

Experimental implementation of the tau adaptive algorithm. This is an approximate method that builds off the tau_leaping method. It self-adapts the value of tau over the course of the simulation. For systems with a small number of molecules, it will be similar in speed to the direct method. For systems with a large number of molecules, it will be much faster than the direct method.

class cayenne.algorithms.tau_adaptive

Runs the adaptive tau leaping simulation algorithm [2]

Parameters:
  • react_stoic ((ns, nr) ndarray) – A 2D array of the stoichiometric coefficients of the reactants. Reactions are columns and species are rows.
  • prod_stoic ((ns, nr) ndarray) – A 2D array of the stoichiometric coefficients of the products. Reactions are columns and species are rows.
  • init_state ((ns,) ndarray) – A 1D array representing the initial state of the system.
  • k_det ((nr,) ndarray) – A 1D array representing the deterministic rate constants of the system.
  • hor – A 1D array of the highest order reaction in which each species appears.
  • nc (int) – The criticality threshold. Reactions with that cannot fire more than nc times are deemed critical.
  • epsilon (float) – The epsilon used in tau-leaping, measure of the bound on relative change in propensity.
  • max_t (float) – The maximum simulation time to run the simulation for.
  • max_iter (int) – The maximum number of iterations to run the simulation for.
  • volume (float) – The volume of the reactor vessel which is important for second and higher order reactions. Defaults to 1 arbitrary units.
  • seed (int) – The seed for the numpy random generator used for the current run of the algorithm.
  • chem_flag (bool) – If True, divide by Na (Avogadro’s constant) while calculating stochastic rate constants. Defaults to False.
Returns:

  • t (ndarray) – Numpy array of the time points.

  • x (ndarray) – Numpy array of the states of the system at times in t.

  • status (int) – Indicates the status of the simulation at exit.

    1 - Succesful completion, terminated when max_iter iterations reached.

    2 - Succesful completion, terminated when max_t crossed.

    3 - Succesful completion, terminated when all species went extinct.

    -1 - Failure, order greater than 3 detected.

    -2 - Failure, propensity zero without extinction.

    -3 - Negative species count encountered

References

[2]Cao, Y., Gillespie, D.T., Petzold, L.R., 2006. Efficient step size selection for the tau-leaping simulation method. J. Chem. Phys. 124, 044109. doi:10.1063/1.2159468