This is just an idea (I won’t be working on it due to other priorities).
igraph has some very basic stochastic SIR modelling functionality.
It would be a cool extension to have a function which runs a similar stochastic model on an arbitrary network, but with an arbitrary number of species (instead of S, I, R) and arbitrary transitions between them, which can also be encoded as a weighted directed network (now we have S --beta--> I and I --gamma--> R). It should allow retrieving individual node states.
Of course, there should be similar (and much more flexible) tools already. Perhaps this? https://github.com/biomathman/pycellerator But I expect they don’t work on a network. They probably consider all-to-all interactions.