You can compute a $LDL^T$ factorization of your covariance matrix, i.e., a factorization in which $L$ is lower triangular with ones on the diagonal and $D$ is diagonal. Then you can interpret every possible "replacement" as a rank-2 update to $D$, so it is simple to compute its determinant. The dominant cost would be solving a linear system with $L$, which costs $O(n^2)$.

Then when your "replacement" take place for real you have to update the LDL^T factorization, and there are algorithms to do that in $O(n^2)$ as well.

I am looking now at the manual page for Matlab; they do not have algorithms for updating LDL but only for Cholesky, but it is essentially the same thing up to diagonal scaling.