I am trying to solve the following system in R, using ode23 (pracma package).

The model I am trying to solve:

http://i.stack.imgur.com/EY0x6.png

I have these initial values and want to run the model for 400 days:

d0 = 0.003, 
d1 = 0.008, 
d2 = 0.05, 
d3 = 1, 
ry = 0.008, 
ay = 1.6/100, 
by = 10/750, 
cy = 100, 
u = 4 * 10^−8, 
y0 = 2.5 * 10^5.

This is my code:

install.packages("deSolve")
library(deSolve)

model <- function(t,x,params){

  y0 <- x[1]
  y1 <- x[2]
  y2 <- x[3]
  y3 <- x[4]

  ry <- params[1]
  mu <- params[2]
  d0 <- params[3]
  ay <- params[4]
  d1 <- params[5]
  by <- params[6]
  d2 <- params[7]
  cy <- params[8]
  d3 <- params[9]

  m <- rep(0,4)
  m[1] = ((ry*(1-mu)) - d0) * y0
  m[2] = (ay * y0) - (d1 * y1)
  m[3] = (by * y1) - (d2 * y2)
  m[4] = (cy * y2) - (d3 * y3)

  return(m)

}

ay = 1.6 / 100
d1 = 0.008
by = 10 / 750
d2 = 0.05
cy = 100
d3 = 1

y_0 = 250000
y_1 = ay*y_0 / d1
y_2 = by*y_1 / d2
y_3 = cy*y_2 / d3

x <- ode23(model, y0 = c(y_0, y_1, y_2, y_3), t0=0,tf=400, params = c(0.008,4*10^-8,0.003,1.6/100,0.008,10/750,0.05,100,1))

plot(x$t,log10(x$y[,1]),xlab='Time (days)', ylab='', type="l",lwd=4)
plot(x$t,log10(x$y[,2]),xlab='Time (days)', ylab='', type="l",lwd=4)
plot(x$t,log10(x$y[,3]),xlab='Time (days)', ylab='', type="l",lwd=4)
plot(x$t,log10(x$y[,4]),xlab='Time (days)', ylab='', type="l",lwd=4)

Equations are from http://ped.fas.harvard.edu/files/ped/files/nature05b_0.pdf?m=1425933222 The plots I get should be similar to the ones in Figure 4.a. in the above link. But it's completely the opposite, all I'm getting are the red lines that are in Figure 4.c. How can I get the plots in Figure 4.a.?