Use ggplot2 with the geom_histogram() function and just set binwidth to whatever you like (this is likely how that figure was made). For fitting, see the R functions glm.nb() and such.

I think the general strategy should be that you find the set of parameters for you chosen model that best fits the data, for example via maximum likelihood. Then for plotting you plug these parameters in your model to get the corresponding values at a grid of points.

In a R these could be done like this:

require(MASS)
set.seed(2)
dat<- rnbinom(1000, mu = 5, size = 3) ## Obs data points
par.nb<- fitdistr(dat, "negative binomial")
par.pois<- fitdistr(dat, "Poisson")

at<- min(dat):max(dat)
p.pois<- dpois(at, par.pois$estimate)
p.nb<- dnbinom(at, size= par.nb$estimate[1], mu= par.nb$estimate[2])

hist(dat, breaks= 25, freq= FALSE, col= 'grey80',
    border= 'white', ylim= c(0, 0.2))
lines(names(table(dat)), table(dat)/length(dat), lwd= 2, col= 'grey30')
lines(at, p.pois, lwd= 2, col= 'blue')
lines(at, p.nb, lwd= 2, col= 'red')
legend('topright', lwd=2, cex= 0.8, col= c('grey30', 'red', 'blue'),
    legend= c('Obs.', 'Fitted neg. binom.', 'Fitted Poisson'))