b = 0.3;
d = 0.2;

pop[0] = 10;
tot[0] = pop[0];
cells = Table[1,{pop[0]}];

tmax = 50;
Do[
pop[t] = pop[t-1];
tot[t] = tot[t-1];
Do[
If[cells[[i]] == 1,
If[RandomReal[]<b,
AppendTo[cells,1];
pop[t] = pop[t]+1;
tot[t] = tot[t]+1
];
If[RandomReal[]<d,
cells[[i]] = 0;
pop[t] = pop[t]-1
]
],
{i,1,tot[t-1]}
],
{t,1,tmax}
];

plotdata = Table[{t,pop[t]},{t,0,tmax}];
plot1 =
ListPlot[
plotdata,
AxesOrigin-> {0,0},
AxesLabel->{"t","pop"},
PlotRange->{{0,51.5},{0,2550}},
PlotStyle->Red,
Filling->Axis,
FillingStyle->Red
]

Epop[t_] := (1-d+b)^t * pop[0];
plot2= 
Plot[Epop[t],{t,0,tmax},PlotStyle->Blue,PlotRange->{{0,51.5},{0,2550}}
];

Show[
plot1,plot2
]

Nruns = 1000;
Do[
pop[0,r] = 10;
tot[0] = pop[0,r];
cells = Table[1,{pop[0,r]}];
Do[
pop[t,r] = pop[t-1,r];
tot[t] = tot[t-1];
Do[
If[cells[[i]] == 1,
If[RandomReal[]<b,
AppendTo[cells,1];
pop[t,r] = pop[t,r]+1;
tot[t] = tot[t]+1
];
If[RandomReal[]<d,
cells[[i]] = 0;
pop[t,r] = pop[t,r]-1
]
],
{i,1,tot[t-1]}
],
{t,1,tmax}
],
{r,1,Nruns}
];

plotdata1 = 
Flatten[
Table[
{t,pop[t,r]},{t,0,tmax},{r,1,Nruns}
],
1
];
Histogram3D[plotdata1]

tallydata =
Tally[plotdata1];
plotdata2 = 
Table[Flatten[tallydata[[i]]],{i,1,Length[tallydata]}];
ListPlot3D[
Select[plotdata2,0<#[[1]]≤tmax &],
PlotRange->All
]

ListPlot3D[
Select[plotdata2,0<#[[1]]≤10&],
PlotRange->All
]

plotdata3 = Table[pop[tmax,r],{r,1,Nruns}];
plotdata3 = Sort[plotdata3];
Print["2.5%ile = ",plotdata3[[Floor[0.025*Nruns]]]];
Print["50%ile = ",plotdata3[[Floor[0.5*Nruns]]],"; expected pop = ",Epop[tmax]];
Print["97.5%ile = ",plotdata3[[Ceiling[0.975*Nruns]]]];
Histogram[plotdata3]

Vpop[t_] := 
Module[{r,beta,N0},
r = 1-d+b;
beta = r - ((1-d)^2+b^2)/2;
N0 = pop[0];
beta*r^(t-1)*N0*(1-r^t)/(1-r)
];

Print["mean -/+ std.dev. = (",Epop[tmax]-Sqrt[Vpop[tmax]],",",Epop[tmax]+Sqrt[Vpop[tmax]],")"];
Print["mean -/+ 2 std.dev. = (",Epop[tmax]-2Sqrt[Vpop[tmax]],",",Epop[tmax]+2Sqrt[Vpop[tmax]],")"];