Suppose there are $n$ variables. If we could easily integrate probabilities over all $n$ variables at once, we would find the exact value for $p$, and would have no need of bounds.

So we assume only two-variable integrals, and then the following algorithm gives non-trivial results. We use $P[S_{k,B,M}:\mu,\Sigma]$ for the probability that $Z_1$ is the max of $M\cdot Z+B$, where \begin{align} B& \text{ is }n\times 1\text{ and constant }\\ M& \text{ is }n\times k\text{ and constant }\\ Z& \text{ is }k\times 1\text{ and random }\\ \mu& \text{ is }k \times 1\text{ and the mean of }Z\\ \Sigma& \text{ is }k\times k\text{ and the covariance of }Z\\ \end{align}

Claim: We will show how, given $k,B,M,\mu,\Sigma$, we can find parameters $k-1,A,L,\lambda,\Pi$ such that (*): $$P[S_{k,B,M}:\mu,\Sigma]\, \le\, P[S_{k-1,A,L}:\lambda,\Pi]$$

Repeating this procedure will eventually give us a two-variable distribution such that: $$P[S_{n,0,1}:\mu,\Sigma]\, \le\, P[S_{2,A,L}:\lambda,\Pi]$$

By hypothesis we can evaluate that right-hand probability to get the desired upper bound. A similar algorithm would also work with $\ge$'s to get the desired lower bound instead. These suffice to answer the question from the claim. $\square$

Proof of Claim:

Let $\sigma_i := \sqrt{\Sigma_{ii}}$ be the standard deviation of $Z_i$.

Let $\theta_{ij} := \arccos(\Sigma_{ij}/\sigma_i\sigma_j) = \arccos(\text{cor}(Z_i,Z_j))$

Let $\ell$ be such that $\mu_\ell$ is the lowest $\mu$ among all $\mu_i$'s except $\mu_1$.

Let $k$ be such that $\mu_k$ is the lowest $\mu$ among all $\mu_i$'s except $\mu_1$ and $\mu_\ell$.

Let $h$ be the number distinct from $k$ and $\ell$ which maximizes $\cos(\theta_{hk}+\theta_{h\ell})$. (For lower bounds, find $h$ which maximizes $\cos(|\theta_{hk}-\theta_{h\ell}|)$ instead.)

Let $\Sigma'$ have all the same entries as $\Sigma$ except that $\Sigma'_{k\ell}=\sigma_k \sigma_\ell \cos(\theta_{hk}+\theta_{h\ell}) \le \Sigma_{k\ell}$

Then $P[S_{n,B,M}:\mu,\Sigma]\le P[S_{n,B,M}:\mu,\Sigma']$, since $\Sigma'$ has a lower correlation.

Furthermore $\Sigma'$ has the lowest possible correlation between $Z_k$ and $Z_\ell$ given their correlations with $Z_h$, so it describes a distribution where $Z_h$, $Z_k$, and $Z_\ell$ are coplanar. Suppose $Z_\ell = a Z_h + b Z_k + c$.

Let $A_i = B_i + cM_{ik}$.

Let $L_i$ (the ith row of the matrix) be the same as $M_i$ but with $L_{ih}=M_{ih}+aM_{i\ell}$, $L_{ik}=M_{ik}+bM_{i\ell}$ and $L_{i\ell}$ omitted.

Let $\lambda$ be $\mu$ but with the entry for $Z_\ell$ omitted.

Let $\Pi$ be $\Sigma'$ but with the row and column for $Z_\ell$ omitted.

This yields the starred claim, as desired, and with that provides the desired algorithm. $\square$

Suppose there are $n$ variables. If we could easily integrate probabilities over all $n$ variables at once, we would find the exact value for $p$, and would have no need of bounds. Instead we will provide bounds by a single two-variable integral. 

We go from $r=n$ to $r=2$, where \begin{align} B^n=0,\text{ and generally }&B^r \text{ is }n\times 1\ \text{ and constant }\\ M^n=1,\text{ and generally }&M^r \text{ is }n\times r\ \text{ and constant }\\ Y^n=Z,\text{ and generally }&Y^r \text{ is }r\times 1\ \text{ and random }\\ \lambda^n=\mu,\text{ and generally }&\lambda^r \text{ is }r \times 1\ \text{ and the mean of }Y,\\ \Pi^r=\Sigma,\text{ and generally }&\Pi^r \text{ is }r\times r\ \text{ and the covariance of }Y, \end{align}

At each step we define $$p=P[Y^r_1=\max M^rY^r+B^r, Y^r \sim N(\lambda^r,\Pi^r)]$$ (we omit superscripts on scalars). We will get upper [lower] bounds by choosing variables so that the $p$'s increase [decrease] as $r$ decreases.

Let $\theta^r_{ij} := \arccos(\Pi^r_{kij}/\sqrt{\Sigma_{ii}\Sigma_{jj}}) = \arccos(\text{cor}(Y^r_i,Y^r_j))$

Let $\ell$ be such that $\lambda^r_\ell=\min_{i>1}\lambda^r_i$.

Let $k$ be such that $\lambda^r_k=\min_{i>1,i\neq \ell}\lambda^r_i$.

Let $h$ be the number distinct from $k$ and $\ell$ which maximizes $\cos(\theta_{hk}+\theta_{h\ell})$. (For lower bounds, find $h$ which maximizes $\cos(|\theta_{hk}-\theta_{h\ell}|)$ instead.)

Let $T$ be the same as $\Sigma$ except that $T_{k\ell}=\sigma_k \sigma_\ell \cos(\theta_{hk}+\theta_{h\ell}) \le \Sigma_{k\ell}$

$T$ has the lowest possible correlation between $Y_k$ and $Y_\ell$ given their correlations with $Y_h$, so it describes a distribution where $Y_h$, $Y_k$, and $Y_\ell$ are linearly dependent.

Let $a,b,c$ be such that $Y_\ell = a Y_h + b Y_k + c$ under $T$.

Let $B^{r-1} = B^r + cM^r_{\cdot k}$.

Let $M^{r-1}_i$ (the ith row of the matrix) be the same as $M^r_i$ but with $M^{r-1}_{ih}=M^r_{ih}+aM^r_{i\ell}$, $M^{r-1}_{ik}=M^r_{ik}+bM^r_{i\ell}$ and $M^{r-1}_{i\ell}$ omitted.

Let $\lambda^{r-1}$ be $\lambda^r$ but with the entry for $Z_\ell$ omitted.

Let $\Pi^{r-1}$ be $T$ but with the row and column for $Z_\ell$ omitted.

Then the new $p$ will be larger than the old $p$, as desired. $\square$