One way to do this will be profile likelihood. If we have a parameter vector $\psi$, profile likelihood is usually calculated for one of the components of $\psi$, but it can be defined for any parametric function of $\psi$. Below is a definition, suppose $L(\psi)$ is the likelihood function and interest (or focus) is on a scalar function $\theta = \theta(\psi)$, then
$$ L_P(\theta) = \max_{\{\psi\colon \theta(\psi)=\theta \}} L(\psi)$$
The implementations of profile likelihood in R (elsewhere?) is not of these generality, so let us make it "by hand".

The model is 
$$ \DeclareMathOperator{\P}{\mathbb{P}} 
p_x= \P(Y=1 \mid X=x)= \frac1{1+e^{-\beta_0 - \beta_1 x}} $$ and the interest parameter $\theta$ is
$$ \theta = p_{0.75} - p_{0.25} $$
It doesn't look promising to try to solve the optimization symbolically, so we try numerically. This is a first attempt, so maybe we can do better. First, a plot of the (negative) profile likelihood for $\theta$, using the data simulated in the question:

[![negative profile loglik for theta][1]][1]

the two blue lines are cutoffs for confidence intervals of 95 and 99%, respectively, based on quantiles from the reference chi-square distribution with 1 df. R code is below:

```r
### First run code from question
library(bbmle)

make_negloglik <- function(y, x) {
   n <- length(y)
   stopifnot( n == length(x) )
   Vectorize( function(beta0, beta1) 
       sum(ifelse(y==0, log1p(exp(beta0  +  beta1*x)),
                                   log1p(exp(-beta0 - beta1*x)))) )
    }

negloglik <- make_negloglik(y, x)

mod.bb <- bbmle::mle2(negloglik,  start=list(beta0=-2, beta1=4))

mod.prof <- bbmle::profile(mod.bb)

plot(mod.prof) # Not shown 

grid <- expand.grid(beta0=seq(-2.8, -0.5, len=100),
                    beta1=seq(1.8, 7.1, len=100))
grid$negloglik <-  with(grid, negloglik(beta0, beta1)) 

P <- function(beta0, beta1, x) 1/( 1  + exp( -beta0 -beta1 * x))

theta <- function(beta0, beta1) P(beta0, beta1, 0.75) - P(beta0, beta1, 0.25)

### Adding theta as a column to data.frame grid:

grid$theta <- with(grid, theta(beta0, beta1))

profile_negloglik <- function(grid) {
    rt <- with(grid,  range(theta))
    seq_theta <- seq(rt[1], rt[2], len=201)
    delta <- diff(seq_theta[1:2])
    npl <- numeric(length=length(seq_theta))
    for (t in seq_along(seq_theta)) {
        tt <- seq_theta[t]
        npl[t] <- with(grid, min(grid[ (tt-delta/2 <= theta) & (theta <= tt + delta/2),
                                      "negloglik" ]))
        }
    return(data.frame(theta=seq_theta, npl=npl))
    }

npl_frame <- profile_negloglik(grid)

npl_min <- with(npl_frame, min(npl))

library(ggplot2)

ggplot(npl_frame, aes(theta, npl))  +  geom_line(color="red")  +
    ggtitle("Profile negative loglikelihood for theta")  +
    geom_hline(yintercept=npl_min)  +
    geom_hline(yintercept=npl_min +  
    qchisq(0.95, 1)/2, color="blue")  +
    geom_hline(yintercept=npl_min +  
    qchisq(0.99, 1)/2, color="blue")   +  ylim(52, 70)   
```  
The idea of the code is:      
* Define a rectangle in parameter space given by individual 99% confidence intervals (calculated by profiling with the R package `bbmle`)      
* use `expand.grid` to cover the rectangle    
* add to the grid data frame a column with the negative loglikelihood, another column with $\theta$       
* Find the range of $\theta$ and subdivide it in many small intervals     
* For each of the intervals, find the minimum negative log likelihood over the interval, and associate that with the midpoint
* finally, plot this as an approximation of the negative profile loglikelihood function of $\theta$. 

As a comparison, let us also calculate an approximate 96% confidence interval using the delta method. Calculations in R: 
```r
 theta_grad <- deriv(expression(  1/( 1  + exp( -beta0 -beta1 * 0.75)) 
              -  1/( 1  + exp( -beta0 -beta1 * 0.25))), 
              c("beta0", "beta1"), function.arg=TRUE)

grad <- theta_grad(coef(model)[1], coef(model)[2])

grad
(Intercept) 
  0.4880566 
attr(,"gradient")
           beta0      beta1
[1,] -0.05555914 0.06582565

grad <- attr(grad, "gradient")
V <- vcov(model)

theta.se <- sqrt( grad %*% V %*% t(grad)  )

 ( CI <- c(0.4881 -2*theta.se, 0.4881  + 2*theta.se ) )
[1] 0.3154351 0.6607649
```       
which is quite close to the profile interval.  

  [1]: https://i.sstatic.net/VWaBi.png