Stochastic Chemical Reaction Simulator

A Python project to simulate stochastic chemical reaction networks using:

• Gillespie Stochastic Simulation Algorithm (SSA)
• Tau-leaping (approximate, faster)

The simulator reads a reaction system from a YAML configuration file and produces
a time series of molecule counts (output.csv) and a plot of the dynamics (output_plot.png), like this:

Features

• Stochastic simulation of reaction networks
• Two algorithms:
  • Exact Gillespie SSA
  • Approximate Tau-leaping with a fixed time step
• Configuration via a simple YAML file:
  • Initial species counts
  • Reactions (reactants, products, rate)
  • External events at specific times (optional)
  • Algorithm choice and random seed
• Exports results to CSV and visualizes them with matplotlib
• Easy to adapt to new models by changing only the configuration file

Tech Stack

• Language: Python 3.10+
• Core libraries: NumPy, pandas, matplotlib, PyYAML
• Paradigm: discrete-event stochastic simulation
• Input format: YAML configuration files

---

Project Structure

.
├─ main.py
├─ base_simulator.py
├─ gillespie_ssa.py
├─ tau_leaping.py
├─ reaction.py
├─ state.py
├─ event.py
├─ data/
│  └─ input.yaml
├─ output.csv          # generated
└─ output_plot.png     # generated

---

Requirements

• Python 3.10+

Python libraries:

• numpy
• pandas
• matplotlib
• pyyaml

Install the dependencies with:

pip install numpy pandas matplotlib pyyaml

(Using a virtual environment is recommended but not required.)

---

How to Run

1. Clone or download this repository.

2. Install the required packages:

   pip install numpy pandas matplotlib pyyaml

3. Make sure the configuration file exists at:

   data/input.yaml
   

4. Run the simulator:

   python main.py

The program will:

• run the simulation until t_end = 0.1
• save the log to output.csv
• generate and save the plot output_plot.png

---

Example Configuration (data/input.yaml)

algorithm: gillespie  # 'gillespie' or 'tau'
events: []
reactions:
  - products:
      Y1: 2
    rate: 10.0
    reactants:
      Y1: 1

  - products:
      Y2: 2
    rate: 0.01
    reactants:
      Y1: 1
      Y2: 1

  - products:
      Z: 1
    reactants:
      Y2: 1
    rate: 1.0

seed: 20
species:
  Y1: 100
  Y2: 100
  Z: 0

Using Tau-leaping

To run the Tau-leaping algorithm instead of Gillespie SSA, change the algorithm field in data/input.yaml:

algorithm: tau

Species

Defines the initial number of molecules for each species (Y1, Y2, Z).

Reactions

Each reaction includes:

• reactants — molecules consumed
• products — molecules produced
• rate — reaction rate constant

The simulator calculates propensities for each reaction based on the current state and decides which reaction occurs next.

Events (optional)

Define external interventions at specific times:

events:
  - time: 0.05
    changes:
      Y1: 10   # add 10 molecules of Y1
      Y2: -5   # remove 5 molecules of Y2

---

Algorithms

Gillespie SSA (Exact)

1. Compute propensities for all reactions.
2. Sample the time to the next reaction from an exponential distribution.
3. Choose which reaction fires based on probabilities.
4. Update molecule counts and repeat.

Pros: Exact stochastic behavior
Cons: Slow for systems with many frequent reactions

Tau-leaping (Approximate)

1. Use a fixed step τ.
2. Estimate how many times each reaction occurs in τ (Poisson distribution).
3. Apply those reactions simultaneously.
4. Advance time by τ.

Pros: Much faster
Cons: Approximate for large τ values

Author

Shahd Amer
Master’s student in Computer Science (Software Programme), University of Pisa (UniPi)

Feel free to open issues, share feedback, or suggest improvements.

Project Context

This project was developed as an individual assignment to practice stochastic simulation algorithms and scientific Python programming.