for func in init_reaction.rate_laws[0].expression.functions:
                for species in func.expression.species:
                    init_species_counts[species.id] = species.distribution_init_concentration.mean
                for obs in func.expression.observables:
                    for species in obs.expression.species: