I think that there is indeed some possibility to lower the bound, and this is something I've looked at seriously a few times. I spent a semester (in 2019) with the Computational Number Theory Group here at BYU trying to do it, without success. Let me outline some of the difficulties we found.

Issue 1. The main technique used in establishing the current bound is optimizing a quotient of integrals. This integral expression was first developed by Maynard and Tao (independently). The optimized integral quotient, in turn, can be approximated very well by finding the largest eigenvalue of a product of certain matrices, $M_1 M_2^{-1}$, related to those integral expressions.

To get better approximations, one needs to increase the size of the matrices. One of the roadblocks in improving knowledge of numerics is that the sizes of these matrices can be too large to store reasonably. I've spent some time optimizing code in Mathematica, storing the matrices as sparse association lists. One still runs into storage issues very easily, even in small dimensions. I'm happy to share this code with anyone who is interested. (Feel free to email me. Disclaimer: There are no comments explaining the code.)

That said, the real bottleneck is the next issue.

Issue 2. Generating of the sparse entries that give rise to the eigenvalue problem takes time. To get to the current bound of $246$ in the Polymath 8b project we needed to slightly generalize the integral quotient mentioned above to an "epsilon-enlarged" region. This adds some complexity to the computations used to form the matrices. This epsilon-enlarged region is still amenable to quick computation, but the bound of $246$ might be the best one can get here.

The polymath project did go beyond the epsilon-enlarged region, using "vanishing marginals", and this led (after a huge amount of effort) to the prime gap bound of $6$ (under a Generalized Elliot-Halberstam conjecture). Unfortunately, these vanishing marginal conditions do not seem to be amenable to quick computation in medium-sized dimensions (or even small dimensions).

In the polymath 8b paper, there are different versions of the epsilon-enlargement idea, and perhaps one of these can immediately give an improved upper bound, but I've not been able to make that work.

To sum up: I believe that there is some small amount of improvement possible to the current bound, and this could easily be done if someone were to improve the epsilon-enlarged region in a way that is amenable to quick computations.

I think that there is indeed some possibility to lower the bound, and this is something I've looked at seriously a few times. I spent a semester (in 2019) with the Computational Number Theory Group here at BYU trying to do it, without success. Let me outline some of the difficulties we found.

Issue 1. The main technique used in establishing the current bound is optimizing a quotient of integrals. This integral expression was first developed by Maynard and Tao (independently). The optimized integral quotient, in turn, can be approximated very well by finding the largest eigenvalue of a product of certain matrices, $M_1 M_2^{-1}$, related to those integral expressions.

To get better approximations, one needs to increase the size of the matrices. One of the roadblocks in improving knowledge of numerics is that the sizes of these matrices can be too large to store reasonably. I've spent some time optimizing code in Mathematica, storing the matrices as sparse association lists. One still runs into storage issues very easily, even in small dimensions. I'm happy to share this code with anyone who is interested. (Feel free to email me. Disclaimer: There are no comments explaining the code.)

That said, the real bottleneck is the next issue.

Issue 2. Generating the sparse entries that give rise to the eigenvalue problem takes time. To get to the current bound of $246$ in the Polymath 8b project we needed to slightly generalize the integral quotient mentioned above to an "epsilon-enlarged" region. This adds some complexity to the computations used to form the matrices. This epsilon-enlarged region is still amenable to quick computation, but the bound of $246$ might be the best one can get here.

The polymath project did go beyond the epsilon-enlarged region, using "vanishing marginals", and this led (after a huge amount of effort) to the prime gap bound of $6$ (under a Generalized Elliot-Halberstam conjecture). Unfortunately, these vanishing marginal conditions do not seem to be amenable to quick computation in medium-sized dimensions (or even small dimensions).

In the polymath 8b paper, there are different versions of the epsilon-enlargement idea, and perhaps one of these can immediately give an improved upper bound, but I've not been able to make that work.

To sum up: I believe that there is some small amount of improvement possible to the current bound, and this could easily be done if someone were to improve the epsilon-enlarged region in a way that is amenable to quick computations.