How to solve a 2D ODE system of the form $\frac{\text{d}\vec{x}}{\text{d}t}=(M-\Delta e^{-\lambda t})\vec{x}+\vec{x}_0$ Let $x=x(t)$ and $y=y(t)$ for $t>0$, with $x(0)=y(0)=1$.
\begin{align}
\frac{\text{d}x}{\text{d}t}
%%%
&=
%%%
-\left( A + \alpha e^{-\lambda t}\right)x + By + x_0\\
%%%
\frac{\text{d}y}{\text{d}t}
%%%
&=
%%%
Cx -\left( D + \beta e^{-\lambda t}\right)y + y_0
\end{align}
where $A,B,C,D,\alpha,\beta,x_0,y_0\in \mathbb{R}_{>0}$. We can rewrite this in the following vector format
\begin{align}
\frac{\text{d}\vec{x}}{\text{d}t}
%%%
&=
%%%
\left(
M
-
\Delta e^{-\lambda t}
\right)
\vec{x}
+
\vec{x}_0
\end{align}
where
\begin{align}
%%%
\vec{x} = 
\left(\begin{array}{c}
x \\
y
\end{array}\right),\qquad
%%%
M = 
\left(\begin{array}{cc}
-A  &  B \\
 C  & -D
\end{array}\right),\qquad
%%%
\Delta = 
\left(\begin{array}{cc}
\alpha & 0 \\
0      & \beta
\end{array}\right),\qquad
%%%
\vec{x}_0 = 
\left(\begin{array}{c}
x_0 \\
y_0
\end{array}\right)
%%%
\end{align}
where $\vec{x}(0)=(1,1)^T$.
The only issue here, that separates it from a standard non-homogeneous linear system of ODEs (see here), is the presence of the exponential vector term.
This 2D system is a reduction of a larger ODE system and this is as simple as I can get it. There does exist a small parameter ($\alpha\simeq\beta\sim O(1/\varepsilon)$, where $0<\varepsilon\ll1$ and all other terms are $O(1)$), however I would rather consider an exact solution as asymptotic analysis on similar systems as this has proved not useful past small/intermediate time in the parameter spaces of interest.
If I make any progress on this I will add it below, if this system has been solved before I would really appreciate a reference, or a reference to a more general solution.
EDIT 1
Ok, so thank you all for your help. I have consulted with my colleagues, discussed their ideas and your suggestions, we believe we have proven that you can't find a solution in terms of elementary functions, even in a reduced case (which is what I need). 
With the aid of some biological arguments I can restrict my parameter regimes of interest to ones that allow the following expression to be surprisingly accurate:
$$y = \frac{1}{1+\gamma e^{-\lambda t}}$$
where $\gamma>0$ and can be expressed in terms of the other coefficients in this model. This reduces the system to a single ODE, as follows:
\begin{align}
\frac{\text{d}x}{\text{d}t}
%%%
&=
%%%
-\left( A + \alpha e^{-\lambda t}\right)x + \frac{B}{1+\gamma e^{-\lambda t}} + x_0
\end{align}
with $x(0)=1$.
NOTE 1:
If we define $\bar{M}$ as follows
\begin{align}
\bar{M} = \left(\begin{array}{cc}
-(A-1)  &  B/\delta \\
\delta C  & -(D-1)
\end{array}\right)
\end{align}
where $\delta\in\mathbb{R}_{>0}$, then the eigenvalues of $\bar{M}$ are $-|\lambda|$ and $-|\mu|$, where $0<\lambda\ll\mu$. If we consider the similar system:
\begin{align}
\frac{\text{d}\vec{z}}{\text{d}t}
%%%
&=
%%%
\bar{M}
\vec{z}
\end{align}
where $\vec{z}=(z_1(t),z_2(t))^T$ and $\vec{z}(0)=(0,1)^T$. Then this system has approximate solution
\begin{align}
\vec{z} &\simeq
\left(\begin{array}{c}
\alpha \\
\beta
\end{array}\right)
e^{-\lambda t}
\end{align}
This solution does not hold for $t\ll1$, as can be seen from the initial conditions, however for all intensive purposes it can be treated as exact.
NOTE 2:
Here is a plot using suitable parameter values:

 A: Solution below is only valid for $\alpha = \beta$, i.e. $\Delta=\alpha I$
A system like this can be solved using a special technique when the matrix $A = M - \Delta e^{\lambda}t$ commutes with it's own integral, i.e. $A(t)\cdot \int_a^tA(\tau)d\tau = \int_a^tA(\tau)d\tau \cdot A(t)$. If you can make this happen ($\alpha = \beta$ for example does it), then you can use that $\Phi(t) = e^{\int_a^tA(\tau)d\tau}$ is a fundamental matrix of solutions for $X^\prime = AX$. You can then use variation of parameters to calculate the solution for the inhomogeneous system. 
See here for some more info on this process: https://www.math24.net/linear-systems-differential-equations-variable-coefficients/
See also this paper for condition on when a matrix commutes with its integral: http://www.jstor.org/stable/2034625

It turns out that $\alpha = \beta$ is a necessary and sufficient condition for $A$ commuting with its integral. If we let $\beta = \alpha$ so that $\Delta = \alpha I$, Then $\Phi = e^\Gamma = \sum_{k=0}^\infty \frac{\Gamma^k}{k!}$, where $\Gamma = \begin{bmatrix} -[At+\frac{\alpha}{\lambda}(e^{\lambda t} - 1)] & B \\ C & -[Dt+\frac{\alpha}{\lambda}(e^{\lambda t} - 1)] \end{bmatrix}$. This gives the homogeneous solution $X_h = e^\Gamma C$ with $C = (c_1,c_2)^T$.
Using variation of parameters, we can get the inhomogeneous solution, by writing $X_i = e^\Gamma K(t)$ and determining $K^\prime (t) = e^{-\Gamma}X_0$, where $X_0 = (x_0,y_0)^T$. This gives the full solution to be $X = X_h + X_i = e^\Gamma(C+\int e^{-\Gamma}X_0)$. Note that $\Gamma =\Gamma(t)$ so the integral is probably not calculable in elementary functions. $C$ can be determined from initial conditions.

With the edit, the problem is reduced to solving a single, linear first order ODE. However, this still does not have a closed form solution. By using an integrating factor (or Mathematica for the lazy) it can be shown that the solution will be
$x(t) =$
$$-e^{\frac{\alpha e^{-\lambda t}}{\lambda}-\frac{\alpha}{\lambda}-A t} \left(e^{\frac{\alpha}{\lambda}} \int_1^0
   \frac{ e^{-\frac{\alpha}{\lambda e^{\lambda \xi } }+A \xi -\lambda \xi } \left(x_0 \gamma e^{\lambda \xi }+x_0 e^{2 \lambda \xi }+B
   e^{2 \lambda \xi }\right)}{\gamma+e^{\lambda \xi }} \, d\xi -e^{\frac{\alpha}{\lambda}} \int_1^t \frac{e^{-\frac{\alpha}{\lambda e^{\lambda \xi
   } }+A \xi -\lambda \xi } \left(x_0 \gamma e^{\lambda \xi }+x_0 e^{2 \lambda \xi }+B e^{2 \lambda \xi }\right)}{\gamma+e^{\lambda \xi
   }} \, d\xi -1\right)$$
For simple parameter choices, this can be easily written in the form of an exponential integral.
A: 
EXPONENTIAL SERIES
Coefficients of the issue system have the imaginary period $\dfrac{2\pi i}\lambda,$ and that allows us to look for the required solution in the form of
$$\vec z = \genfrac{(}{)}{0}{0}{x}{y} = \sum\limits_{n=0}^\infty{\vec{c_n} e^{-\lambda nt}},\tag{A}$$
then
$$\dfrac{d\vec z}{dt} = \genfrac{(}{)}{0}{0}{x}{y} = -\lambda\,\sum\limits_{n=1}^\infty{n\vec{c_n} e^{-\lambda nt}}.$$
Since
$$\dfrac {d\vec z}{dt} = \left(M - \Delta e^{-\lambda t}\right)\vec z + \vec z_0,$$
where 
$$\vec {z_0} = \genfrac{(}{)}{0}{}{x_0}{y_0},$$
one can get:
$$\left(M + \gamma\,nE\right)\vec{c_{n}} = \Delta\vec{c_{n-1}},\ n\in\mathbb N,$$
$$\vec{c_n} = \left(\prod\limits_{i=1}^n \left((M + \gamma\,nE)^{-1} \Delta\right)\right)\vec{c_0},$$
$$(M - \Delta)\sum\limits_{n=0}^\infty \vec{c_n}= - \vec{z_0}.$$
So
$$\boxed{\vec{c_0} = - \left(E + \sum\limits_{n = 1}^\infty P_n\right)^{-1}(M - \Delta)^{-1}z_0,\quad \vec{c_n} = P_n \vec{c_0},\,}\tag{B}$$
where
$$P_n = \prod\limits_{i=1}^n \left((M + \gamma\,nE)^{-1}\Delta\right),\quad n\in\mathbb N.\tag{C}$$

LOOKING FOR THE COMMON SOLUTION
At first, let us using the nice idea of $\,\mathbf{user121049}$.
Differentiating of the first equation leads to the system
\begin{cases}
\dfrac{d^2x}{dt^2} = \alpha\lambda e^{-\lambda t}\, x - \left(A + \alpha e^{-\lambda t}\right)\dfrac{dx}{dt} + B\dfrac{dy}{dt} + x_0\\
\dfrac{dy}{dt} = Cx - \left(D + \beta e^{-\lambda t}\right)y + y_0\\
By = \dfrac{dx}{dt} + \left(A + \alpha e^{-\lambda t}\right)x - x_0\\
\end{cases}
$$\dfrac{d^2x}{dt^2} = \alpha\lambda e^{-\lambda t}\, x - \left(A + \alpha e^{-\lambda t}\right)\dfrac{dx}{dt} + B\left(Cx + y_0\right) - \left(D + \beta e^{-\lambda t}\right)\left(\dfrac{dx}{dt} + \left(A + \alpha e^{-\lambda t}\right)x - x_0\right) + x_0,
$$
$$\dfrac{d^2x}{dt^2} + 2(G + \gamma e^{-\lambda t})\dfrac{dx}{dt} + H(t)x  =  {By_0  - \left(D + \beta e^{-\lambda t} + 1\right)x_0},\tag1
$$
where
$$G = \frac{A+D}{2},\quad \gamma = \frac{\alpha + \beta}{2},\quad H(t) = AD - BC + (A\beta + (D - \lambda)\alpha)e^{-\lambda t} +\alpha\beta e^{-2\lambda t}.\tag2$$
Let
$$x = uv,\quad x' = u'v + uv'',\quad x''= u''v +2u'v' + uv'',$$
then LSH of $(1)$ looked up
$$u''v + 2u'v' + uv'' + 2(G+ \gamma e^{-\lambda t}) (u'v + uv')  H(t)uv.$$
Zeroing of the $u'$ factor can be achieved when
$$v' + \left(G + \gamma e^{-\lambda t}\right) v = 0,\quad v(t) = \exp\left(-Gt + \frac\gamma\lambda e^{-\lambda t}\right),\tag3$$
and then
$$v'(t) = -\left(G + \gamma e^{-\lambda t}\right)v,$$
$$v''(t) = \left(\gamma\lambda e^{-\lambda t} + \left(G + \gamma e^{-\lambda t}\right)^2\right)v,$$
$$x = uv,\tag4$$
$$x' = \left(u'- \left(G + \gamma e^{-\lambda t}\right)u\right)v,$$
$$x'' = \left(u'' - 2\left(G + \gamma e^{-\lambda t}\right)u' + \left(\gamma\lambda e^{-\lambda t} + \left(G + \gamma e^{-\lambda t}\right)^2\right)u\right)v,$$
and, taking in account $(2),$ equation $(1)$ can be reconstructed to the form of
$$\dfrac {d^2u}{dt^2} + \left(\gamma\lambda e^{-\lambda t} - \left(G + \gamma e^{-\lambda t}\right)^2 - H(t)\right)\,u = {\left(By_0  - \left(D + \beta e^{-\lambda t} + 1\right)x_0\right)e^{Gt}\exp\left(-\frac\gamma\lambda e^{-\lambda t}\right)},$$
$$\dfrac {d^2u}{dt^2} - \left(-\frac12\lambda(\alpha + \beta)e^{-\lambda t} + \frac14\left(A + D + (\alpha + \beta) e^{-\lambda t}\right)^2\right)\,u - \left(AD - BC + (A\beta + (D - \lambda)\alpha)e^{-\lambda t} +\alpha\beta e^{-2\lambda t}\right)\,u = {\left(By_0  - \left(D + \beta e^{-\lambda t} + 1\right)x_0\right)e^{Gt}\exp\left(-\frac\gamma\lambda e^{-\lambda t}\right)},$$
$$\dfrac {d^2u}{dt^2} - \left(BC + \frac14\left(A-D + (\alpha - \beta)e^{-\lambda t}\right)^2\right)\,u = {\left(By_0  - \left(D + \beta e^{-\lambda t} + 1\right)x_0\right)\exp\left(\frac{A + D}2t\right)\exp\left(-\frac{\alpha + \beta}{2\lambda} e^{-\lambda t}\right)},\tag5$$
The homogenius analogue of $(5)$ is known as Hill equation, and
its solutions can be expressed throw the special functions (Whittaker function). 
On the other hand, substitution 
$$u'(t) = u(t)w(t)\tag6$$
transforms the Hill equation to Riccati equation in the form of
$$\dfrac {dw}{dt} + w^2 = BC + \frac14\left(A-D + (\alpha - \beta)e^{-\lambda t}\right)^2,\tag7$$
and all my further attempts also have not bring serious result.
These circumstances make applying of the constant variation method too hard.
