Fitting SIR model with 2019-nCoV data doesn't conververge I am trying to calculate the basic reproduction number $R_0$ of the new 2019-nCoV virus by fitting a SIR model to the current data. My code is based on https://arxiv.org/pdf/1605.01931.pdf, p. 11ff:
library(deSolve)
library(RColorBrewer)

#https://en.wikipedia.org/wiki/Timeline_of_the_2019%E2%80%9320_Wuhan_coronavirus_outbreak#Cases_Chronology_in_Mainland_China
Infected <- c(45, 62, 121, 198, 291, 440, 571, 830, 1287, 1975, 2744, 4515)
day <- 0:(length(Infected)-1)
N <- 1400000000 #pop of china
init <- c(S = N-1, I = 1, R = 0)
plot(day, Infected)


SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, { dS <- -beta * S * I
  dI <- beta * S * I - gamma * I
  dR <- gamma * I
  list(c(dS, dI, dR))
  })
}

RSS.SIR <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = day, func = SIR, parms = parameters)
  fit <- out[ , 3]
  RSS <- sum((Infected - fit)^2)
  return(RSS)
}

lower = c(0, 0)
upper = c(0.1, 0.5)

set.seed(12)
Opt <- optim(c(0.001, 0.4), RSS.SIR, method = "L-BFGS-B", lower = lower, upper = upper)
Opt$message
## [1] "NEW_X"

Opt_par <- Opt$par
names(Opt_par) <- c("beta", "gamma")
Opt_par
##      beta     gamma 
## 0.0000000 0.4438188

t <- seq(0, 100, length = 100)
fit <- data.frame(ode(y = init, times = t, func = SIR, parms = Opt_par))
col <- brewer.pal(4, "GnBu")[-1]
matplot(fit$time, fit[ , 2:4], type = "l", xlab = "Day", ylab = "Number of subjects", lwd = 2, lty = 1, col = col)
points(day, Infected)
legend("right", c("Susceptibles", "Infecteds", "Recovereds"), lty = 1, lwd = 2, col = col, inset = 0.05)


R0 <- N * Opt_par[1] / Opt_par[2]
names(R0) <- "R0"
R0
## R0 
##  0

I also tried fitting with GAs (as in the paper), also to no avail.
My question
Am I making any mistakes or is there just not enough data yet? Or is the SIR model too simple? I would appreciate suggestions on how to change the code so that I get some sensible numbers out of it.
Addendum
I wrote a blog post based on the final model and current data: Epidemiology: How contagious is Novel Coronavirus (2019-nCoV)?
 A: You might be experiencing numerical issues due to the very large population size $N$, which will force the estimate of $\beta$ to be very close to zero. You could re-parameterise the model as
\begin{align}
{\mathrm d S \over \mathrm d t} &= -\beta {S I / N}\\[1.5ex]
{\mathrm d I \over \mathrm d t} &= \beta {S I / N} - \gamma I \\[1.5ex]
{\mathrm d R \over \mathrm d t} &= \gamma I \\
\end{align}
This will make the estimate of $\beta$ larger so hopefully you'll get something more sensible out of the optimisation.
In this context the SIR model is useful but it only gives a very crude fit to these data (it assumes that the whole population of China mixes homogenously). It's perhaps not too bad as a first attempt at analysis. Ideally you would want some kind of spatial or network model that would better reflect the true contact structure in the population. For example, a metapopulation model as described in Program 7.2 and the accompanying book (Modeling Infectious Diseases in Humans and Animals, Keeling & Rohani). However this approach would require much more work and also some data on the population structure. An approximate alternative could be to replace the $I$ in $\beta SI/N$ (in both of the first two equations) with $I^\delta$ where $\delta$, which is probably $<1$, is a third parameter to be estimated. Such a model tries to capture the fact that the force of infection on a susceptible increases less than linearly with the number of infecteds $I$, while avoiding specification of an explicit population structure. For more details on this approach, see e.g. Hochberg, Non-linear transmission rates and the dynamics of infectious disease, Journal of Theoretical Biology 153:301-321.
A: Because the population of china is so huge, the parameters will be very small.
Since we are in the early days of the infection, and because N is so big, then $S(t)I(t)/N \ll 1$.  It could me more reasonable to assume that at this stage of the infection, the number of infected people is approximately exponential, and fit a much simpler model.
A: This is only marginally related to the detailed coding discussion, but seems highly relevant to the original question concerning modeling of the current 2019-nCoV epidemic.     Please see arxiv:2002.00418v1 (paper at https://arxiv.org/pdf/2002.00418v1.pdf ) for a delayed diff equation system ~5 component model, with parameter estimation and predictions using dde23 in MatLab.   These are compared to daily published reports of confirmed cases, number cured, etc.  To me, it is quite worthy of discussion, refinement, and updating.    It concludes that there is a bifurcation in the solution space dependent upon the efficacy of isolation, thus explaining the strong public health measures recently taken, which have a fair chance of success so far.
A: what do you think about putting the initial number of infectious as an addition parameter in the optimization problem otherwise the fitting need to start with the initial condition.
