Edit: I have recently created two Jupyter Notebooks that illustrate this behaviour and let you play around with some actual matrices and actual signals.
I find understanding MDCT easiest if we define our transforms as matrix operations.
DFT
In case of a DFT such a matrix would look like
$$
F_M = \left( \sqrt{\frac{1}{M}} e^{-2 \pi i k n / M} \right)_{k,n=0 \dots M-1}
\in \mathbb{C}^{M \times M}
$$
$F_M \approx$ 
and our transform would become a matrix multiplication
$$
X = F_M x
\qquad{} x,X \in \mathbb{C}^M
$$
where $X$ and $x$ are column vectors of size $M$. The inverse of the matrix $F$ exists, so synthesis is possible
$$
x = F^{-1}_M X
\qquad{} x,X \in \mathbb{C}^M
$$
which means (note we are using the unitary matrix definition here)
$$
F_M F^{-1}_M = I
$$
Wide DFT
If we now just defined a "wide DFT" matrix to be larger in one dimension, lets say $M \times 2M$.
$$
G_{2M} = \left( \sqrt{\frac{1}{M}} e^{-2 \pi i k n / M} \right)_{k=0 \dots M-1 \quad{} n=0 \dots 2M-1}
\in \mathbb{C}^{M \times 2M}
$$
this is equivalent to having longer basis functions.
$G_{2M} \approx$ 
This means for $2M$ samples $x$ we get $M$ coefficients $\hat{X}$.
$$
\hat{X} = G_{2M} x
\qquad{} x \in \mathbb{C}^{2M}
\quad{} X \in \mathbb{C}^M
$$
Doing that, you might also notice is that there are now repititions in our matrix as our basis function $e^{-2 \pi i k n / M}$ starts to repeat after $M$ samples. It looks like two DFT matrices stacked next to each other, which we can represent as a folding matrix followed by DFT
$$
G_{2M} = F_M
\begin{bmatrix}
I_M & I_M
\end{bmatrix}
$$
$G_{2M} \approx$ $=$ $\times$ 
Where $I$ is the identity matrix. Obviously an inverse of $\begin{bmatrix}I_M & I_M\end{bmatrix}$ doesn't exist, so an inverse of $G$ doesn't exist either.
If we now look at the complete transform
$$
\hat{X} = F_M
\begin{bmatrix}
I_M & I_M
\end{bmatrix} x
\qquad{} x \in \mathbb{C}^{2M}
\quad{} X \in \mathbb{C}^M
$$
we see that always the sum of exactly two samples of $x$ are being transformed. This is called time domain aliasing.
We can do this for any size and the folding matrix will just repeat, e.g.
$$
G_{4M} = F_M
\begin{bmatrix}
I_M & I_M & I_M & I_M
\end{bmatrix}
$$
$G_{4M} \approx$ $=$ $\times$ 
and we will end up with the sum of more and more samples being transformed together, producing more and more aliasing.
DCT-IV
We can define DCT-IV as a matrix, too.
$$
C_M = \left( {\sqrt{\frac{2}{M}} \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{1}{2}\right)\right]} \right)_{k,n = 0 \dots M-1}
\in \mathbb{R}^{M \times M}
$$
$C_M \approx$ 
Interestingly, the matrix is orthogonal and its own inverse
$$
C_M = C^{-1}_M = C^T_M
$$
Wide DCT-IV
Again, we can define a wide DCT-IV of size $M \times 2M$ by just extending the basis functions.
$D_{2M} \approx$
However, this time the repititions in the basisfunction are different, and thus the folding matrix is:
$$
D_{2M} = C_M
\begin{bmatrix}
I_M & -J_M
\end{bmatrix}
\in \mathbb{R}^{M \times 2M}
$$
$D_{2M} \approx$ = $\times$ 
where $J_M$ is the $M \times M$ antidiagonal matrix. If we repeat this further, we get
$$
D_{4M} = C_M
\begin{bmatrix}
I_M & -J_M & -I_M & J_M
\end{bmatrix}
\in \mathbb{R}^{M \times 4M}
$$
$D_{4M} \approx$ 
= $\times$ 
after which our folding matrix will start to repeat. You will discover this is exactly the repetition pattern of the image on Wikipedia
We will need this pattern, so lets define it
$$
Q =
\begin{bmatrix}
I_M & -J_M & -I_M & J_M
\end{bmatrix}
\in \mathbb{R}^{M \times 4M}
$$
MDCT
Analysis
$$
X_k = {\sqrt{\frac{2}{M}} \sum^{2M-1}_{n=0} h(n) x_n \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{M}{2}+\frac{1}{2}\right)\right]} \qquad{} k = 0 \dots M-1
$$
If we look at the definition of the MDCT we note that we are transforming $n = 0 \dots 2M$ samples to $k = 0 \dots M$ coefficients, so we are using a "wide matrix" of size $M \times 2M$.
$$
E_{2M} = \left( {\sqrt{\frac{2}{M}} \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{M}{2}+\frac{1}{2}\right)\right]} \right)_{k = 0 \dots M-1 \quad{} n=0 \dots 2M-1}
\in \mathbb{R}^{M \times 2M}
$$
$$
X_k = E_{2M} x
\qquad{} x \in \mathbb{R}^{2M}
\quad{} X \in \mathbb{R}^M
$$
But also, if we compare the basis function to DCT-IV, we notice they are almost the same, but with an offset of $M/2$ for $n$. This means we are using a "shifted view" of DCT-IV matrix, shifted to the right:
$$
E_{2M} = \left( {\sqrt{\frac{2}{M}} \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{1}{2}\right)\right]} \right)_{k = 0 \dots M-1 \quad{} n=\frac{M}{2} \dots \frac{M}{2}+2M-1}
\in \mathbb{R}^{M \times 2M}
$$
We can express such a shifted view as a cyclic permutation matrix $P\in \mathbb{R}^{4M \times 4M}$, but only if our basis function was cyclic, too. We know $D_{4M}$ is expanded to the point where it becomes cyclic, so we can express $E_{2M}$ using $D_{4M}$ and using an extraction matrix and a permutation matrix:
$$
E_{2M} =
D_{4M}
P
\begin{bmatrix}
I_{2M} \\
0
\end{bmatrix}
$$
$E_{2M} \approx$ 
$=$ $\times$ 
$\times$ 
but since we know how to represent wide DCT-IV by a regular DCT-IV like
$$
D_{4M} = C_M Q
$$
we can write
$$
E_{2M} = C_M Q P
\begin{bmatrix}
I_{2M} \\
0
\end{bmatrix}
$$
$E_{2M} \approx$ $=$ $\times$ $\times$ $\times$
We don't need all of those shift matrices individually, so let's multiply them together into one matrix $R$:
$$
R = Q P
\begin{bmatrix}
I_{2M} \\
0
\end{bmatrix}
$$
$E_{2M} \approx$ $=$ $\times$ 
And end up with
$$
X = C_M R x
$$
This means MDCT is a permutation and superposition and then a DCT-IV of our input samples.
Synthesis
For reconstruction we need to form a linear combination of basis vectors, which can be written as
$$
\hat{x} = (C_M R)^T X = R^T C_M^T X
$$
And substituing our forward transform for $X$ we get
$$
\hat{x} = R^T C_M^T C_M R x
$$
we know
$$
C_M^TC_M = I_M
$$
but
$$
R^T R \neq I_M
$$
instead
$R^T R \approx$ 
But if we now overlap+add two frames with just the right hopsize
we get
perfect reconstruction in the overlapping part.
This means, if we first transformed our signal in windows with half a window overlap, all the aliasing produced in this overlap will exactly cancel out, making the entire transform perfectly reconstructing.
Caveats/Further reading
Windowing
All of this is obviously disregarding windowing. But if you look at the permutation matrix and look at what sample is superpositioned with what other sample it should be at least somewhat intuitively visible where the Princen Bradley condition comes from.
Nicer Permutation
Some authors prefer the MDCT to be defined as
$$
X_k = {\sqrt{\frac{2}{M}} \sum^{2M-1}_{n=0} h(n) x_n \cos\left[\frac{\pi}{M} \left(k+\frac{1}{2}\right) \left(n+\frac{M}{2}-\frac{1}{2}\right)\right]} \qquad{} k = 0 \dots M-1
$$
(note the minus). If you then solve all of the remaining math, you will find a permutation matrix that looks like
$\hat{R} \approx$ 
as opposed to
$R \approx$ .
This doesn't change any of the remaining explanations (TDAC still works etc.), but you can see that the samples in the center of your frame are now 'nicer': they're positive valued and transformed forwards, as opposed to negative valued and transformed backwards.
Math for overlap add
If you feel that last bit became a bit too "hand wavey", that's because it is. The traditional math of transforms or matrices starts to break down if you start talking about overlapping transforms. A very neat way of dealing with them are Polyphase Matrices in the z-Domain. They're hard to wrap your head around at first though.
