Can the infinite sum $\sum_{n=0}^\infty {2^n \sum_{k=0}^n (-1)^k \frac{ {{n}\choose{k}}}{ (n+k)! }}$ be simplified? 
EDIT2:
   I have made this related question for general $n$ (in this question $n=2$).

Can the sum
$$S = \sum_{n=0}^\infty  {2^n \sum_{k=0}^n (-1)^k \frac{ {{n}\choose{k}}}{ (n+k)! }}$$
be simplified?
I think the inner sum $ b_n = \sum_{k=0}^n (-1)^k \frac{ {{n}\choose{k}}}{ (n+k)! }$ equals the integral
$$\int_{0}^1 (1-x_1) \int_{0}^{1-x_1} (1-x_2) \dots \int_{0}^{1-\sum_{j=1}^{n-1} x_j} (1-x_n) dx_n\dots dx_2dx_1$$
since I calculated it for some values of $n$ and found OEIS A006902. Is this equality true?
EDIT:
Hypothesis: $S$ equals the limit (as $M\to \infty$) of the sum of the first column of the matrix
$$
A_M = -M^2
\begin{pmatrix}
-M^2 &  &  \\
2M-1 & -M^2  &  \\
2M-3 & 2M-1 & -M^2  \\
\vdots & \vdots & \vdots & \ddots  &  \\
5 & 7 & 9 & \dots & -M^2  & \\
3 & 5 & 7 & \dots & 2M-1 & -M^2
\end{pmatrix}^{-1}
$$
The matrix $A_M$ has $-M^2$ on the diagonal (or equivalently $1$ if the multiplier $-\frac{1}{M^2}$ is put inside the inverse and multiplied into the matrix) and then on the lower diagonals the odd numbers starting from $2M-1$ down to $3$.
I guess (from calculating these for some values of $M$) that the first column of $A_M$ is $\left(1, \frac{\alpha_2}{M^2}, \frac{\alpha_3}{M^4},\dots \frac{\alpha_M}{M^{2(M-1)}} \right)$, where the $\alpha$'s are some natural numbers (for example for $M=5$ they are $(1,9, 256, 7004, 182836)$ ($\alpha_1 = 1$ always)). A WA calculation. It also seems that  $\alpha_2 = 2M-1$ and $\frac{\alpha_j}{\alpha_{j-1}} \to M^2$.
 A: Let $\ell^\infty = \left\{ (x_n)_{n\ge 0} : \sup_n |x_n| < \infty \right\}$ be the space of bounded sequences over $\mathbb{C}$. 
Equipped with sup-norm $\|x\|_\infty = \sup_n|x_n|$, it is a Banach space. Let $B(\ell^\infty)$ be the collection of bounded linear operators on $\ell^\infty$. Equipped with operator norm, it is a Banach algebra. We will use these facts to justify the algebraic operations below.
Let $\epsilon \in \ell^\infty$ be the sequence $\left( \frac{1}{n!} \right)_{n\ge 0}$.
Define a linear operator $L$ over $\ell^\infty$ by shifting the entries of a sequence to the left. More precisely,
$$\ell^\infty \ni x = (x_0,x_1,\ldots) \quad\mapsto\quad Lx = (x_1,x_2,\ldots) \in \ell^\infty$$
It is easy to see $L \in B(\ell^\infty)$ with operator norm $\|L\| = 1$.
For any $x \in ( 0, 4 )$, let $a = \sqrt{\frac1x - \frac14} > 0$
and $\mu = x \left(\frac12 + a i\right)$.
Consider following sequence:
$$s(x) = (s_0(x),s_1(x),\ldots)\quad\text{ where }\quad
s_\ell(x) = \sum_{n=0}^\infty x^n \left[\sum_{k=0}^n (-1)^k \binom{n}{k}\frac{1}{(n+k+\ell)!}\right]$$
In terms of $L$ and $\epsilon$, we can formally express $s(x)$ as
$$s(x) 
= \sum_{n=0}^\infty x^n \left[\sum_{k=0}^n (-1)^k \binom{n}{k} L^k \right] L^n \epsilon
= \sum_{n=0}^\infty (xL(1-L))^n \epsilon\tag{*1}
$$
When $|x| < \frac12$, we have
$$\| xL(1-L) \| \le |x|\|L\| + |x|\|L\|^2 \le 2|x| < 1
$$
and the expansion on RHS of $(*1)$ converges. This means for $x \in (0,\frac12)$, $s(x)$ equals to
$$\begin{align}
s(x) = (s_0(x),s_1(x),\ldots)   
&= \frac{1}{1 - xL(1-L)} \epsilon
= \frac1x \frac{1}{\left(L - \frac12\right)^2 + a^2} \epsilon\\
&= \frac{i}{2ax}\left[\frac{1}{L - \frac12 + a i} - \frac{1}{L - \frac12 - ai}\right]
\epsilon\\
&= \frac{1}{2axi}\left[\frac{\mu}{1 - \mu L} - \frac{\bar{\mu}}{1 - \bar{\mu}L}\right]
\end{align}
\tag{*2}
$$
Please note that when $x \in (0,\frac12)$, $|\mu| = \sqrt{x}< 1$ and the two inverses on last line are well defined.
For any $f = (f_0,f_1,\ldots) \in \ell^\infty$, let $f(z)$ be the power series $\sum_{\ell=0}^\infty f_\ell z^\ell$. 
Since $\|f\|_\infty < \infty$, $f(z)$ converges for $|z| < 1$. 
In particular, we have $\epsilon(z) = \sum_{\ell=0}^\infty \frac{z^\ell}{\ell!} = e^z$.
Let $g = \frac{1}{1 - \mu L}\epsilon = \sum_{k=0}^\infty \mu^{k}L^k \epsilon$. 
The corresponding function $g(z)$ equals to
$$g(z) = \sum_{k=0}^\infty\sum_{\ell=0}^\infty
\frac{\mu^k z^\ell}{(k+\ell)!}
= 
\sum_{n=0}^\infty \frac{1}{n!}\left(
\frac{z^{n+1} - \mu^{n+1}}{z - \mu}\right)
= \frac{ze^z -\mu e^{\mu}}{z-\mu}
$$
By comparing the coefficients of $z^0$ on both sides, we obtain 
$$\frac{1}{1-\mu L}\epsilon = (e^{\mu},\ldots)$$
By a similar argument, we have
$$\frac{1}{1-\bar{\mu} L}\epsilon = (e^{\bar{\mu}},\ldots)$$
Substitute these into $(*2)$, we find for $x \in (0,\frac12)$
$$\begin{align}s_0(x)
&= \frac{1}{2axi}\left(\mu e^{\mu} - \bar{\mu}e^{\bar{\mu}}\right)
= e^{\frac{x}{2}}\left[\cos(xa) + \frac{\sin(xa)}{2a}\right]
\\
&=
e^{\frac{x}{2}}\left[
\cos\left(\frac12\sqrt{x(4-x)}\right)
+ \sqrt{\frac{x}{4-x}}\sin\left(\frac12\sqrt{x(4-x)}\right)
\right]
\end{align}
$$
The RHS of above expression defines a function analytic near the origin. Since the singularity nearest to the origin is located at $x = 4$, we can analytic continue above expression to all $x \in \mathbb{C}$ with $|x| < 4$. 
In particular, at $x = 2$, this reduces to the sum at hand:
$$\sum_{n=0}^\infty 2^n \left[\sum_{k=0}^n (-1)^k \binom{n}{k}\frac{1}{(n+k)!}\right] = s_0(2) = e(\cos(1) +\sin (1))
\approx 3.756049227094727
$$
