In the end, I came up with the following approach, which is basically a hierarchical model.

The historical data was used as "true" distribution of successful catches in 7 tasks. By sampling from this distribution, I create a set of the true abilities of n mice. Assuming that each mouse has a random component in its coordinative task and won't always perform according to its true ability, I used each mouse's true ability as proportion p in a binomial distribution to sample successes in k=7 trials. Due to the treatment, each mouse's individual p is then lowered by 25% (of course, I could have assumed a random component on the treatment effect, too, but I thought that could be a slight overkill). Finally, a Wilcoxon signed-rank test is used to compare the number of successes before and after treatment.

library(coin) # since this package can handle ties in the data
y <- c(rep(0,2),rep(1,3),rep(2,9),rep(3,17),rep(4,19),rep(5,30),rep(6,15),rep(7,5)) # historical data

nsim=1e4
n=14
power=0
p.value<-c()

while(power<0.8){
  n=n+1
  for(i in 1:nsim){
    smp <- sample(y, n, replace=TRUE)
    before.treat <- rbinom(n=n,size=7,prob=smp/7)
    after.treat <- rbinom(n=n,size=7,prob=(smp/7)*0.7)
    if(all((before.treat-after.treat)==0)){
      p.value[i] <- 1 # if there are no differences at all, we want a non-significant p-value
    } else { 
      test <- wilcoxsign_test(before.treat ~ after.treat)
      p.value[i] <- pvalue(test)
    }
  }
  power<-sum(p.value<0.05)/nsim
  print(paste("N=",n,", power=",power))
}

> [1] "N= 15 , power= 0.7182"
> [1] "N= 16 , power= 0.7413"
> [1] "N= 17 , power= 0.7829"
> [1] "N= 18 , power= 0.8077"

I would be happy to take any commments on this approach. Particularly, I noted that the Monte Carlo simulations seem very sensitive to the random seed used and larger simulations do take a lot of time, so I am a bit worried about the validity of the results.