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 intents and purposes it can be treated as exact.
NOTE 2:
Here is a plot using suitable parameter values:
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.