We introduce expandable partial propensity direct method (EPDM) - a new exact stochastic simulation algorithm suitable for systems involving many interacting molecular species. The algorithm is especially efficient for sparsely populated systems, where the number of species that may potentially be generated is much greater than the number of species actually present in the system at any given time. The number of operations per reaction scales linearly with the number of species, but only those which have one or more molecules. To achieve this kind of performance we are employing a data structure which allows to add and remove species and their interactions on the fly. When a new specie is added, its interactions with every other specie are generated dynamically by a set of user-defined rules. By removing the records involving the species with zero molecules, we keep the number of species as low as possible. This enables simulations of systems for which listing all species is not practical. The algorithm is based on partial propensities direct method (PDM) by Ramaswamy et al. for sampling trajectories of the chemical master equation.