### Metapopulation model function
metapopmodel = function(alpha.gamma, beta.gamma, 
                        alpha.phi, beta.phi1, beta.phi2, beta.phi3, beta.phi4, 
                        alpha.disp, beta.non.urban, beta.urban 
                        ){
  # Some base parameters
  NSITES = nsites # all sites in the original set (167)
  NYEARS = sim.years # number of years to run the simulation forward
  # Preallocate arrays
  occ = consim = logit.gam = gamsim = logit.phi = phisim = poccsim = water = array(dim=c(NYEARS, NSITES))
  pr.disp = array(dim=c(NYEARS, NSITES, NSITES))
  conwsim = array(dim=c(NYEARS, NSITES, NSITES))
  nocc = vector(length=NYEARS) # to monitor the number occupied each year of the simulation
  # Start the model (parameters at first time-step)
  occ[1, ] = occupancy
  nocc[1] = sum(occ[1, ])
  water[1, ] = rbinom(NSITES, 1, pr.water[, 1]) # presence or absence of water at each site in first year (as observed)
  # Dynamics
  random.year = sample(2:21, 50, replace=T) # random year for sampling water presence or absence
  for(t in 2:NYEARS){
    water[t, ] = rbinom(NSITES, 1, pr.water[, random.year[t]]) # presence or absence of water at each site in each year a random sample from all years
    # Probability of dispersal between all sites, based on distance between them (but limit neighbours to just those within 2 km)
    for(i in 1:NSITES){
      for(j in 1:NSITES){
        pr.disp[t, i, j] = ifelse(dist[i, j] > 2, 0, plogis(alpha.disp + beta.non.urban * non.urban.dist.use[i, j] + 
                                                                      beta.urban * urban.dist.use[i, j]))
      }
    }
    diag(pr.disp[t, , ]) = 0 # set the diagonal to zero, to ensure no self connectivity
    # Calculate contribution of each neighbour to connectivity at each time-step (probability of dispersal multipied by neighbour occupancy)
    conwsim[t, , ] = pr.disp[t, , ] * occ[t-1, ]
    # Calculate connectivity as the cumulative probability of dispersal from all occupied neighbours
    consim[t, ] = 1 - apply(1 - conwsim[t, , ], 2, prod)
    # Specify logistic model for the probability of colonisation 
    logit.gam[t, ] = alpha.gamma + beta.gamma * consim[t, ] 
    gamsim[t, ] = sapply(logit.gam[t, ], plogis) * water[t, ] * present.use # multiply by the binary 'water' variable to set colonisation to zero when no water and set to zero when site not present
    # Specify logistic model for the probability of persistence
    logit.phi[t, ] = alpha.phi + beta.phi1 * chytrid.use + beta.phi2 * type + beta.phi3 * aqveg.use + beta.phi4 * consim[t, ] 
    phisim[t, ] = sapply(logit.phi[t, ], plogis) * water[t, ] * present.use
    # Probability of occupancy in each year a function of occupancy status, persistence and colonisation
    poccsim[t, ] = occ[t-1, ] * phisim[t, ] + (1 - occ[t-1, ]) * gamsim[t, ]
    # Occupancy in each year a Bernoulli variable, with probability poccsim
    occ[t, ] = rbinom(NSITES, 1, poccsim[t, ])
    # Number of sites occupied in each year
    nocc[t] = sum(occ[t, ])
  }
  return(nocc)
}