model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dtA <- rA*A*(1 - (A+alpha*C)/kA)
    dtC <- rC*C*(1 - (C+beta*A)/kC)
    return(list(c(dtA,dtC)))  
  }) 
} 

pA <- c(rA=1.11, kA=104.73)
pC <- c(rC=0.916, kC=60.277) 
p <- c(pA,pC,m=1,alpha=0,beta=0)

aurelia1 <- read.table("data/aurelia1.txt", header=TRUE)
aurelia2 <- read.table("data/aurelia2.txt", header=TRUE)
caudatum1 <- read.table("data/caudatum1.txt", header=TRUE)
caudatum2 <- read.table("data/caudatum2.txt", header=TRUE)
A0 <- c(aurelia1[1,2],0,aurelia2[1,2],aurelia2[1,2])
C0 <- c(0,caudatum1[1,2],caudatum2[1,2],caudatum2[1,2])
s <- c(A=A0[1], C=C0[2]);s

# Check previous fits in absence of competition (alpha=beta=0)
run(18,legend=FALSE)
timePlot(aurelia1,draw=points,add=TRUE,legend=FALSE)
timePlot(caudatum1,draw=points,add=TRUE,colMap=2,legend=FALSE)

# Next use these estimates to fit the two data sets from the mixed experiment:

timePlot(aurelia2,draw=points,legend=FALSE)
timePlot(caudatum2,draw=points,add=TRUE,colMap=2,legend=FALSE)
p["alpha"] <- 1 ; p["beta"] <- 1
free <- c("alpha","beta")
f1 <- fit(list(aurelia2,caudatum2),free=free,add=TRUE,legend=FALSE)
summary(f1)

# Since all parameters are known we can sketch nullclines

p[free] <- f1$par; p
plane(xmax=110,ymax=110)
run(20,0.1,traject=TRUE)

# Finally fit all data together by fixing the initial conditions.
# Aurelia and Caudatum are absent from the first and second expt, respectively:

fixed <- list(A=A0,C=C0)
free <- c("rA","kA","rC","kC","alpha","beta")
f2 <- fit(list(aurelia1,caudatum1,aurelia2,caudatum2),free=free,fixed=fixed,add=TRUE,legend=FALSE)
summary(f2) # Now all parameters are identifiable 

p[free] <- f2$par
f3 <- fit(list(aurelia1,caudatum1,aurelia2,caudatum2),free=free,fixed=fixed,add=TRUE,legend=FALSE,bootstrap=100)
pairs(f3$bootstrap,main="Parameter correlations")