At the start of your derivation you multiply out the brackets $\sum_i (x_i - \bar{x})(y_i - \bar{y})$, in the process expanding both $y_i$ and $\bar{y}$. The former depends on the sum variable $i$, whereas the latter doesn't. If you leave $\bar{y}$ as is, the derivation is a lot simpler, because
\begin{align}
\sum_i (x_i - \bar{x})\bar{y}
&= \bar{y}\sum_i (x_i - \bar{x})\\
&= \bar{y}\left(\left(\sum_i x_i\right) - n\bar{x}\right)\\
&= \bar{y}\left(n\bar{x} - n\bar{x}\right)\\
&= 0
\end{align}
Hence
\begin{align}
\sum_i (x_i - \bar{x})(y_i - \bar{y})
&= \sum_i (x_i - \bar{x})y_i - \sum_i (x_i - \bar{x})\bar{y}\\
&= \sum_i (x_i - \bar{x})y_i\\
&= \sum_i (x_i - \bar{x})(\beta_0 + \beta_1x_i + u_i )\\
\end{align}
and
\begin{align}
\text{Var}(\hat{\beta_1})
& = \text{Var} \left(\frac{\sum_i (x_i - \bar{x})(y_i - \bar{y})}{\sum_i (x_i - \bar{x})^2} \right) \\
&= \text{Var} \left(\frac{\sum_i (x_i - \bar{x})(\beta_0 + \beta_1x_i + u_i )}{\sum_i (x_i - \bar{x})^2} \right), \;\;\;\text{substituting in the above} \\
&= \text{Var} \left(\frac{\sum_i (x_i - \bar{x})u_i}{\sum_i (x_i - \bar{x})^2} \right), \;\;\;\text{noting only $u_i$ is a random variable} \\
&= \frac{\sum_i (x_i - \bar{x})^2\text{Var}(u_i)}{\left(\sum_i (x_i - \bar{x})^2\right)^2} , \;\;\;\text{independence of } u_i \text{ and, Var}(kX)=k^2\text{Var}(X) \\
&= \frac{\sigma^2}{\sum_i (x_i - \bar{x})^2} \\
\end{align}
which is the result you want.
As a side note, I spent a long time trying to find an error in your derivation. In the end I decided that discretion was the better part of valour and it was best to try the simpler approach. However for the record I wasn't sure that this step was justified
$$\begin{align}
& =.
\frac{1}{(\sum_i (x_i - \bar{x})^2)^2}
E\left[\left( \sum_i(x_i - \bar{x})(u_i - \sum_j \frac{u_j}{n})\right)^2 \right] \\
& =
\frac{1}{(\sum_i (x_i - \bar{x})^2)^2} E\left[\sum_i(x_i - \bar{x})^2(u_i - \sum_j \frac{u_j}{n})^2 \right]\;\;\;\;\text{ , since } u_i \text{ 's are iid} \\
\end{align}$$
because it misses out the cross terms due to $\sum_j \frac{u_j}{n}$.
Comparison of $I$ with correlation coefficients is good, but it has its limits. This answer uncovers what those limits are. It derives a tight upper bound for $|I|$ in terms of the weights matrix $W$ and shows that in ordinary applications, where $W$ is symmetric and row-normalized, this bound is $1$ (or less). For the extraordinary weights matrix in the question (which has some zero rows, making it impossible to row-normalize), it shows that $X=\pm(3,-1,-1,-1)$ achieves the most extreme value of $|I|$ possible, which is $3$.
Analysis: Simplifying the Formula for $I$
The repeated appearances of $x_i - \bar x$ in the formula, along with the sum of their squares in the denominator, invite us to re-express $I$ in terms of the standardized variates
$$z_i = \frac{(x_i - \bar x)}{\sqrt{\sum_{i=1}^n (x_i - \bar x)^2}};\quad \bar x = \frac{1}{n}\sum_{i=1}^n x_i,$$
which can be assembled into the column vector $z = (z_1, \ldots, z_n)^\prime$. Because the variables have been standardized, $z$ is a unit vector
$$|z|^2 = z_1^2 + \cdots + z_n^2 = 1.$$
(This differs slightly from the usual statistical concept of "standardization," which recenters and rescales $x$ to have length $\sqrt{n}$.)
Further simplicity is achieved by renormalizing the weights $\omega_{ij}$ so they sum to unity:
$$\sum_{i,j} \omega_{ij} = 1.$$
When $S_0 = \sum_{i,j}\omega_{ij} \ne 0$, this can always be done by dividing every weight by $S_0$. (When the sum of the weights is zero, the formula for $I$ is undefined in the first place.) The renormalized weights form a matrix $W = (\omega_{ij} / S_0) = (w_{ij})$.
In terms of $z$ and $W$,
$$\eqalign{
I &= \sum_{i,j} \frac{n \omega_{ij}}{\sum_{i,j} \omega_{ij}} \left(\frac{(x_i-\bar x)}{\sqrt{\sum_i(x_i-\bar x)^2}}\right)\left(\frac{(x_j-\bar x)}{\sqrt{\sum_i(x_i-\bar x)^2}}\right) \\
&= \sum_{i,j} n w_{ij}z_i z_j = z^\prime (n W z).
}$$
This is identical to the formula for the Pearson correlation
$$\rho(z, w) =z^\prime w$$
between any two standardized variables $z$ and $w$.
Comparison to Correlation Coefficients
The difference between $I$ and a correlation is that $w = nWz$ might not be standardized: $|nWz| = \sqrt{(nWz)^\prime (nWz)}$ is not necessarily equal to $1$. To standardize, we must first subtract the mean of $w = nWz$,
$$\bar{w} = \frac{1}{n}\mathbf{1}^\prime n W z = \mathbf{1}^\prime W z,$$
from each component of $n W z$. Let $V$ be the matrix for this transformation,
$$V:z \to Vz = nWz - \mathbf{1}\bar w = \left(n\mathbb{I}_n - \mathbf{1}\mathbf{1}^\prime\right) W z.$$
(The vector $\mathbf{1}$ is the column vector $(1,1,\ldots,1)^\prime$, so the matrix $\mathbf{1}\mathbf{1}^\prime$ is the $n\times n$ matrix all of whose entries are ones.)
By standardizing $Wz$ in this manner--subtracting $\bar w$ and dividing by the length of the resulting centered vector--we get an honest correlation coefficient whose size cannot exceed $1$ (by virtue of the Cauchy-Schwarz Inequality):
$$1 \ge |\rho(z, Wz)| = \frac{|z^\prime Vz|}{|Vz|}.$$
This can be related to $I$ because $\mathbf{1}^\prime z = 0$, whence $z^\prime V z = z^\prime W z$. Clearing the denominator in the preceding inequality gives
$$|I| = |z^\prime (nW z)| = n|z^\prime V z| \le |V z|.$$
Equality will hold if and only if $Vz \ne 0$ and $Vz$ is parallel to $z$: that is, when $z$ is an eigenvector of $V$ with nonzero eigenvalue. When such an eigenvector exists (and it almost always will), this shows that
The size of $I$ is bounded by the largest eigenvalue (in absolute value) of $V = (n\mathbb{I} - \mathbf{1}\mathbf{1}^\prime)W$ and it can attain this bound.
Whenever we can find a basis in which $V$ is diagonal the largest absolute entry of $V$ is its largest eigenvalue.
In the question the unit-sum-normalized version of the weights matrix is
$$W = \pmatrix{1&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0}$$
whence
$$V = \pmatrix{ 3&0&0&0 \\ -1&0&0&0 \\ -1&0&0&0 \\ -1&0&0&0 }.$$
It is simple to compute that the largest value attained by $|Vz|$ among all centered unit vectors $z$ (that is, $\mathbf{1}^\prime z = 0$ and $|z|=1$) is $3$, uniquely attained at $z=\pm(3,-1,-1,-1)/\sqrt{12}$. Consequently, all we can hope in general for this particular weights matrix $W$ is that $|I| \le 3$.
Resolution of the Original Question
Ordinarily, $W$ has nonnegative entries that originally are row standardized to sum to unity. This automatically makes the sum of each row in the unit-sum-normalized version of $W$ equal to $1/n$. Consequently $n W \mathbf{1} = \mathbf{1}$, making $\mathbf{1}$ an eigenvector of $n W$ with eigenvalue $1$.
It is also well-known that in this case if $W$ is symmetric, all other eigenvectors $z$ will be orthogonal to $\mathbf{1}$--that is, $z$ will be zero-centered--and will have eigenvalue $\lambda$ between $0$ and $1$. Considering any such unit $z$,
$$z^\prime V z = z^\prime\left(n\mathbb{I}_n - \mathbf{1}\mathbf{1}^\prime\right)W z = z^\prime (nW z) = z^\prime (\lambda z) = \lambda z^\prime z = \lambda \le 1.$$
In this ordinary application of Moran's I we may conclude that $\max{|I|} = 1$, as expected. But when $W$ is not symmetric and row-standardized--or cannot possibly be row-standardized because one or more rows are entirely zero--then the maximum has to be computed using the general result attained previously.
Best Answer
It looks like the null model that the expectation is taken over makes this easily proved. The null model is that one picks a random permutation $\pi$ from the set of all permutations uniformly at random and the null model is
$$\mathbb{E}_{\pi} [I] = \frac{N}{w} \frac {\sum_i \sum_j W_{ij}(({\pi x})_i-\bar {\pi x}) (({\pi x})_j-\bar {\pi x})} {\sum_i (({\pi x})_i-\bar {\pi x})^2}$$
We have that $$\bar {\pi x} = \bar {x},$$ $$\sum_i (({\pi x})_i-\bar {\pi x})^2 = \sum_i (x_i-\bar x)^2$$
Further, $$\begin{eqnarray} \mathbb{E}_{\pi}[(({\pi x})_i-\bar {\pi x}) (({\pi x})_j-\bar {\pi x})] &= \mathbb{E}_{\pi}[(({\pi x})_i-\bar {x}) (({\pi x})_j-\bar {x})],\\ &= \frac{1}{N(N-1)} \sum_{i \ne j} (x_i - \bar {x})(x_j - \bar {x}),\\ &= \frac{1}{N(N-1)} [ (\sum_r x_r - \bar {x})^2 - \sum_r (x_r - \bar {x})^2 ],\\ &= \frac{-\sum_r (x_r - \bar {x})^2}{N(N-1)} \end{eqnarray}$$
Thus $$\begin{eqnarray} \mathbb{E}_{\pi} [I] &= \frac{N}{w} \frac {(\sum_i \sum_j W_{ij})(\frac{-\sum_r (x_r - \bar {x})^2}{N(N-1)})} {\sum_i (x_i-\bar x)^2},\\ &= \frac{-1}{N-1}. \end{eqnarray}$$