Since the terms in the sum are decaying very fast, we can restrict the sum to small values of $k$ and use an asymptotic estimate that has a small error for large $n$ and small $k$. In more detail, if
$$a_k = \binom n k k^{k-1} (n-k)^{2n-k},$$
then
$$\frac {a_{k+1}} {a_k} =
\left( 1 + \frac 1 k \right)^{k-1}
\left( 1 - \frac 1 {n-k} \right)^{2n-2k}
\left( 1 - \frac 1 {n-k} \right)^{k-1} < \frac 1 e, \\
\sum_{k = \ln^2 n + 1}^n a_k < a_1 \sum_{k = \ln^2 n}^{n - 1} e^{-k} =
O \left( n^{2n} e^{-\ln^2 n} \right),$$
thus the tail will be negligible compared to all the terms in the expansion $n^{2n} \left( c_0 + \frac {c_1} n + \dots \right)$. On the other hand, for large $n$,
$$a_k = \frac {k^{k-1} e^{-2 k}} {k!} n^{2 n}
\left( 1 - \frac {k^2 - k} {2 n} +
O \left(\frac {k^4} {n^2} \right) \right), \\
\sum_{k=1}^{\ln^2 n} a_k =
n^{2 n} \left(
\sum_{k=1}^{\ln^2 n} \frac {k^{k-1} e^{-2 k}} {k!}
\left( 1 - \frac {k^2 - k} {2 n} \right) +
O \left(\frac {\ln^{10} n} {n^2} \right) \right),$$
once we prove that the $O$ term is uniform in $k$.
Using the same reasoning, we can now extend the summation range on the right-hand side to infinity, as the contribution from the tail is $O(e^{-\ln^2 n})$. This yields the Lambert W-function. The second term contains the factors $k^p$, which correspond to applying the operators $\left( z \frac d {dz} \right)^p$ to the W-function. After simplifications,
$$\sum_{k=1}^n \binom n k k^{k - 1} (n - k)^{2 n - k} \sim
n^{2 n} \left(
\alpha - \frac {\alpha^2 (2 - \alpha)} {2(1 - \alpha)^3 n} + \dots \right), \\
\alpha = -W \!\left(-e^{-2} \right).$$
