Gillespie algorithm

gillespie.py

import collections
import csv
import random

Reaction = collections.namedtuple('Reaction', ('reaction_rate', 'reaction'))


def get_probability_of_reaction(reaction: Reaction, substrates: dict) -> float:
    probability = reaction.reaction_rate
    for element_name, element_value in reaction.reaction.items():
        if element_value < 0:
            number_of_elements_in_simulation = substrates[element_name]
            how_many_elements_needed = -element_value
            probability *= how_many_elements_needed * number_of_elements_in_simulation
    return probability


def chose_next_reaction(reactions: [Reaction], substrates: dict) -> Reaction:
    weights = [get_probability_of_reaction(reaction, substrates) for reaction in reactions]
    return random.choices(reactions, weights=weights)[0]


def update_simulation(substrates: dict, reaction: Reaction):
    for element_name, element_value in reaction.reaction.items():
        substrates[element_name] += element_value
        assert substrates[element_name] >= 0


def step(reactions: [Reaction], substrates: dict):
    reaction = chose_next_reaction(reactions, substrates)
    update_simulation(substrates, reaction)


def simulate(substrates: dict, reactions: [Reaction], steps: int):
    with open('simulation.csv', 'w') as f:
        writer = csv.DictWriter(f, fieldnames=substrates.keys())
        writer.writeheader()
        for _ in range(steps):
            step(reactions, substrates)
            writer.writerow(substrates)


if __name__ == '__main__':
    elements = {
        'A': 10,
        'B': 10,
        'C': 10,
        'D': 10,
        'E': 10,
    }

    reactions = [
        Reaction(reaction_rate=0.1, reaction={'A': -1, 'B': -1, 'C': +1}),
        Reaction(reaction_rate=0.1, reaction={'D': -1, 'A': +1, 'B': +1}),
        Reaction(reaction_rate=0.1, reaction={'E': +1}),
        Reaction(reaction_rate=0.1, reaction={'E': -1, 'C': -1, 'D': +1}),
        Reaction(reaction_rate=0.1, reaction={'A': -1}),
        Reaction(reaction_rate=0.1, reaction={'B': -1}),
        Reaction(reaction_rate=0.1, reaction={'C': -1}),
        Reaction(reaction_rate=0.1, reaction={'D': -1}),
        Reaction(reaction_rate=0.1, reaction={'E': -1}),
    ]
    simulate(elements, reactions, 100)