I would like to simulate birth process with time delay. Here is a reference ref on page 3. Note there is two species whereas in my case I have one species. So only equation (1) and (2) should be considered. My code does not give desired result. Any suggestion.
Here is the code for Birth Death Process without delay(propensity rate depends on current population).
$ \emptyset\xrightarrow[]{\text{$ \lambda$ X}} X\quad \quad$ birth
$ X\xrightarrow[]{\text{$ \mu$ X}} \emptyset \quad \quad$ death
SeedRandom@2; With[{\[Lambda] = 4, \[Mu] = 1, initialPop = 10}, sim = NestList[( a1 = \[Lambda] #[[2]]; a2 = \[Mu] #[[2]]; a0 = Total@{a1, a2}; reaction = 1/a0 Accumulate@{a1, a2}; pos = First@FirstPosition[reaction, _?(RandomReal[] < # &)]; \[CapitalDelta]t = RandomVariate@ExponentialDistribution[a1 + a2]; Which[ pos == 1, {#[[1]] + \[CapitalDelta]t, #[[2]] + 1}, pos == 2, {#[[1]] + \[CapitalDelta]t, #[[2]]  1} ] ) &, {0, initialPop}, 20]]; sim; ListLinePlot[sim, Epilog > {Red, PointSize[Medium], Point[sim]}, Frame > True, PlotTheme > "Detailed", FrameLabel > {"Time", "Population"}, ImageSize > Large, InterpolationOrder > 0]
Now Assume there is delay on birth reaction.
$ \emptyset\xrightarrow[]{\text{$ \lambda$ X}} X\quad \quad$ birth has a delay
$ X\xrightarrow[]{\text{$ \mu$ X}} \emptyset \quad \quad$ death has no delay
 Sample reaction time $ t_1\sim Exp(\lambda X+\mu X)$
 Choose a reaction.
 If it is a death set $ X=X1$ and new time is $ t_1$ .
 If it is a birth sample $ t_d\sim Gamma(4,2)$
 Put $ t_1+t_d$ into queue.
 Sample new reaction time $ t_w\sim Exp(\lambda X+\mu X)$
 If $ t_1+t_d<t_1+t_w$ set $ X=X+1$ . New time $ t_1+t_d$
 If $ t_1+t_d>=t_1+t_w$
 Choose a reaction.
 If it is a death set $ X=X1$ .
 If it is a birth set $ X=X$ . No birth. New time $ t_1+t_w$

Repeat
SeedRandom@2;
With[{[Lambda] = 30, [Mu] = 1, initialPop = 10}, sim = NestList[( a1 = [Lambda] #[3]; a2 = [Mu] #[3]; a0 = Total@{a1, a2}; reaction = 1/a0 Accumulate@{a1, a2}; pos = First@FirstPosition[reaction, _?(RandomReal[] < # &)]; pos2 = First@FirstPosition[reaction, _?(RandomReal[] < # &)]; [CapitalDelta]t = RandomVariate@ExponentialDistribution[a1 + a2]; td = RandomVariate@GammaDistribution[4, 2]; tw = RandomVariate@ExponentialDistribution[a1 + a2];
Which[ pos == 2, {#[[1]] + \[CapitalDelta]t, #[[2]]  1}, pos == 1, Which[ \[CapitalDelta]t + td < \[CapitalDelta]t + tw, {#[[1]] + \[CapitalDelta]t + td, #[[2]] + 1}, \[CapitalDelta]t + td >= \[CapitalDelta]t + tw, Which[ pos2 == 2, {#[[1]] + \[CapitalDelta]t + tw, #[[2]]  1}, pos2 == 1, {#[[1]] + \[CapitalDelta]t + tw, #[[2]]} ] ] ] ) &, {0, initialPop}, 100]];
sim;
ListLinePlot[sim, Epilog > {Red, PointSize[Medium], Point[sim]}, Frame > True, PlotTheme > “Detailed”, FrameLabel > {“Time”, “Population”}, ImageSize > Large, InterpolationOrder > 0]