If you are open to programs outside of R and Python, you could use CAFE 5 (Mendes et al. 2020). Its aim is to infer rates, but it has an option to simulate, as described in the slow start guide:

Simulate - Given specified values, generate an artificial list of gene families that matches the values. To generate a simulation, pass the --simulate or -s parameter. Either pass a count of families to be simulated with the parameter (--simulate=1000) or pass a --rootdist (-f) parameter with a file containing the distribution to match (see [rootdist] for the file format).

P.S: I am confused by what models you request: you say "under a BM or OU model", and then "birth-and-death". I think there is a misunderstanding because the former models are suited for continuous traits while the latter for counts.