The lancet article
The method from the article in the Lancet is very straightforward. They model the spread of the virus by simulating the transmission on the level of each single infected individual. For each infected person they compute (randomly) how many other people they will infect and how long it will take for those others to infect others (or how likely it is that the others get contained due to policy measures like quarantaine).
They use this model to estimate the variation of the potential number of new cases after some given time, by computing the random model a thousand times for each particular set of model parameters. If the number is large then they consider the outbreak uncontrolled and then this stochastic model can be used to express the probability to control the outbreak for a given set of parameters.
So the difference between a deterministic model and stochastic model is the following:
deterministic The virus spreads with a constant number and and speed. For instance each person passes on the virus to two others in a given specific time interval and then the growth will be like 1,2,4,8,16,etc.
stochastic The spread of the virus is random. By how many it increases is random and not every time the same factor. Some people spread a lot others only a little. For instance sometimes a person passes it on to three other people, and sometimes only one (but on average the same namely two). And then the growth will be random it may be high (when it triples a lot) it may be low (when only one person gets it). This randomness is expressed by repeating the model several times and then see how it ends up in all those cases.
The serial interval distribution is explained in Figure 2 a . I haven't read the article in detail, but after a quick scan it seems to me that the serial time is the time between the moments that a person infects another person. The serial time distribution is the distribution of those times. It is not the name for a specific distribution.
Note that the code of the Lancet article is available online. https://github.com/cmmid/ringbp/tree/master/R
The imperial college article
The infections will not continue exponentially. This is only in the very beginning. The reason why the infections rate decreases is because you can't infect a person that is already infected before. So the probability to spread the virus becomes less in time. (also, the infection rate depends on the weather/season as well, sometimes called respiratory season, that is I believe not incorporated in those models)
A well known model that takes this decrease of infection rate into account is the SIR model (and this will already produce graphs like the ones you are looking for). However that model assumes homogeneous mixing which is not very realistic. So the model that they use at imperial college uses many smaller compartments which are schools, workplaces, households and probably some more. Then the spread occurs differently at different levels/distances. You can not infect your housemates when they are already sick, so often it is only one person that infects all the others members in the family (and that person has a relatively high transmission) and the others will spread relatively less (but they can pass on the virus in other places, like school work church etc. If these are not saturated yet.
This is not easy to replicate. What you need to do is model the spatial structure realistically. Like age distribution in households and the networks of who is going to which work/school/church etc. A description of this work is given in one of the references. The model is normally used for influenza. https://www.pnas.org/content/suppl/2008/02/28/0706849105.DC1
I have made a toy model that sort of showcases these effects (but not with realistic distributions). What you get is not exponential growth but instead more something like a power law. The distribution grows in space and spreads at the edge of the infected population. It is a bit like the growth of area of a circle as function of it's circumference. $dA/dt = \text{constant} \times \text{circumference}$ but then for a fractal dimension structure.
The result of the toy model is a curve that is in the beginning exponential (homogeneous mixing growth), but then changes into a power law relationship (growth at the edge of some geometric figure). In any case the growth is not exponential with a continuous rate, but instead the dynamics is changing (in this example the growth is only exponential for the first 5 generations).
# create 500x500 people in matrix
set.seed(1)
L <- 5*10^2
people <- matrix(rep(0,(L)^2),L)
# trackers for the locations of the people that got sick:
# we start with index patient in the middle
orderx <- c(round(L/2))
ordery <- c(round(L/2))
generation <- c(1)
spread <- 0
R0 <- 3
R1 <- 0.25 # a probabiliy to spread the virus on long distance, e.g. due to travel.
##### run the virus ######
# compute probability density function
# for probabilty of spreading out to nearby locations
Lr <- 7
Lspread <- 1+Lr*2
# targets will be in a cube of LrxLr around the patient
targets <- matrix(1:Lspread^2,Lspread)
xt <- matrix(rep(c(1:Lspread)-(Lspread+1)/2,Lspread),Lspread)
yt <- t(xt)
# ps is some probability to get infected as function of distance
ps <- c(exp(-c(Lr:1)*0.2),0,exp(-c(1:Lr)*0.2))
probs <- ps[xt+(Lspread+1)/2]*ps[yt+(Lspread+1)/2]
### plot for visualization of the spread
plot(orderx,ordery,xlim=c(1,L),ylim=c(1,L),
xlab = "", ylab= "",
col=1,bg = 1,cex=0.2,pch=21)
# itterate all the patients untill all have been dealt with
# during this loop the number of patients increases
while (spread < length(generation)) {
spread <- spread + 1
x <- orderx[spread]
y <- ordery[spread]
g <- generation[spread]
# selecting Rn people in the neighbourhood of the patient
# Rn is sampled from a Poisson distribution with mean R0
Rn <- rpois(1,R0)
if (Rn>0) {
sel <- sample(targets,Rn, prob = probs)
xt[sel]
yt[sel]
## this loop picks out the R0 people
## these are gonna become new patients if they are susceptible
for (i in 1:Rn) {
#the modulo is to patch left with right and top with bottom
xq <- (x+xt[sel[i]]-1)%%L+1
yq <- (y+yt[sel[i]]-1)%%L+1
# if the 'target' is not sick yet then add it as new patient
if (people[xq,yq] == 0) {
generation <- c(generation,g+1)
orderx <- c(orderx,xq)
ordery <- c(ordery,yq)
people[xq,yq] <- g+1
colv <- (g+1)/30-floor((g+1)/30)
points(xq,yq,
col=hsv(colv,1,1),bg = hsv(colv,1,1),cex=0.1,pch=21)
}
}
}
### additionally make R1 random people from far away sick
nfar <- rpois(1,R1)
ifar <- 0
while (ifar<nfar) {
ifar = ifar +1
xq <- sample(1:L,1)
yq <- sample(1:L,1)
if ((people[xq,yq] == 0)*(rbinom(1,1,0.1)==1)) {
generation <- c(generation,g+1)
orderx <- c(orderx,xq)
ordery <- c(ordery,yq)
people[xq,yq] <- g+1
colv <- (g+1)/30-floor((g+1)/30)
points(xq,yq,
col=hsv(colv,1,1),bg = hsv(colv,1,1),cex=0.1,pch=21)
}
}
}
# ratio of people that got sick
spread/L^2
# plot the spread in colours
colv <- (generation+1)/40-floor((generation+1)/40)
plot(orderx,ordery,xlim=c(1,L),ylim=c(1,L),
xlab = "", ylab= "",
col=hsv(colv,1,1),bg = hsv(colv,1,1),cex=0.1,pch=21)
# plot the epidemiological curve
I <- sapply(1:50, FUN = function(x) sum(generation == x))
plot(I, log = 'xy',
xlab = "x, generation", ylab = "number of infectious people", type = "l",
ylim = c(1,5*10^4), xlim = c(1,70))
gen <- 1:50
colv <- (gen+1)/40-floor((gen+1)/40)
points(I,pch=21,col = 1, bg = hsv(colv,1,1))
lines((R0+R1)^c(0:50), lty=2)
sm <- 4:50
lines(sm,0.5*sm^3.5, lty = 3)
lines(sm,0.002*sm^6, lty = 4)
legend(1,5*10^4, c(expression((R[0]+R[1])^x),expression(0.5*x^3.5),
expression(0.002*x^6)), lty = c(2,3,4),
xjust = 0, cex = 0.7)
So my colleagues could play around with the parameters and create different scenarios. 
