Let $x=[x_1, x_2, ..., x_{n+1}]$ is a random vector following a multivariate normal distribution;
$$x \sim N(x; \mu, \Sigma)$$
where $\mu$ is a mean vector of x, and $\Sigma$ is a covariance matrix. Only $x_1$ is correlated with the other variables, and the others are independent of each other. Thus, $\Sigma$ is expressed as:
$$\Sigma = \begin{bmatrix}d_1 & c_1 & c_2 & \cdots & c_n \\c_1 & d_2 & 0 & \cdots & 0 \\c_2 & 0 & \ddots& & 0 \\ \vdots & \vdots & & &\vdots \\c_n & 0 & \cdots & & d_{n+1}\end{bmatrix}$$
In this situation, I want to compute a probability that x is in a hyper-rectangle;
$$\int_{\Omega} f(x) dx$$
where $\Omega=\{{ x \mid a_i \leq x_i \leq b_i }\}$
How can I compute that probability? Or, how can we compute an upper or lower bound for that probability?
Let me change a bit your notations: we assume that the covariance of $(X,Y_1,\ldots,Y_n)$ is, by blocks $$\left[\begin{array}{cc}a&c^T\\c&D\end{array}\right]$$ where $c^T=(c_1,\ldots,c_n)$ and $D=\mathrm{diag}(d_1,\ldots,d_n).$ For simplicity I assume that $E(X)=E(Y_i)=0.$ The theory says that $X$ conditioned by $Y$ is Gaussian with mean $$c^TD^{-1}Y=\sum_{i=1}^n\frac{c_i}{d_i}Y_i=m(Y)$$ with variance independent of $Y$ and equal to $\sigma^2=a-\sum_{i=1}^n\frac{c^2_i}{d_i}.$ Note that $m(Y)$ is Gaussian. As a consequence $$\Pr(x<X<x'|Y)=\Phi(\frac{x'-m(Y)}{\sigma})-\Phi(\frac{x'-m(Y)}{\sigma})=f(Y).$$ Now if $B=\prod_{i=1}^n[y_i,y'_i]$ and if $g(y)=Ce^{-(\sum_{i=1}^n\frac{y_i^2}{2d_i})}$ is the density of Y, then $$\Pr(x<X<x', Y\in B)=\int_Bf(y)g(y)dy.$$