I have the following network SEIR model:
\begin{align} \frac{d}{dt} S_k(t) &= -\lambda kS_k(t) \Theta(t)\\\\ \frac{d}{dt} E_k(t) &= \lambda kS_k(t) \Theta(t) - \varepsilon E_k(t)\\\\ \frac{d}{dt} I_k(t) &= \varepsilon E_k(t) - \gamma I_k(t)\\\\ \frac{d}{dt} R_k(t) &= \gamma I_k(t). \end{align}
where $S_k(t), E_k(t), I_k(t), R_k(t)$ are densities of the population who are susceptible, exposed, infected, and removed at time $t$ for vertices of degree $k$; i.e. $S_k+ E_k+ I_k + R_k=1$. $\Theta(t)$ is the probability of having an infected neighbor given by
$$ \Theta(t) = \sum_{k} \frac{k P(k) I_{k}(t)}{\langle k \rangle}.$$
I have already solved the system but it doesn't really have a closed form solution so all I can do is solve it the best I can and interpret analytically. My problem is, I cannot seem to solve one equation (particularly the 2nd) and was wondering if anyone can help point out what is wrong or give direction to some other approach I could use to solve it.
Here are my solutions thus far:
\begin{align*} \int \frac{d}{dt} S_k(t) &= \int -\lambda kS_k(t) \Theta(t)\\\\ \int \frac{dS_k(t)}{S_k(t)} &= -\lambda k \int \Theta(t) dt\\\\ \ln S_k(t) &= -\lambda k \phi(t)\\\\ S_k(t) &= e^{-\lambda k \phi(t)} \end{align*}
where $\phi$ is the auxiliary function given by $$\phi(t) = \int \Theta(t') dt' = \frac{1}{\langle k \rangle} \int \sum_{k'} kP(k')I_{k'}(t) dt = \frac{\sum_{k'} kP(k')R_{k'}(t)}{\gamma \langle k \rangle}$$ Solution for $R_k(t)$ \begin{align*} \int \frac{d}{dt}R_k(t) &= \int \gamma I_k(t)\\\\ \int dR_k(t) &= \gamma \int I_k(t) dt\\\\ R_k(t) &= \gamma \int I_k(t)dt \qquad \longrightarrow \int I_k(t)dt = \frac{1}{\gamma} R_k(t) \end{align*} Solution for $I_k(t)$ \begin{align*} \int \frac{d}{dt}I_k(t) &= \int \varepsilon E_k(t) - \gamma I_k(t)\\\\ \int dI_k(t) &= \int (\varepsilon E_k(t) - \gamma I_k(t))dt\\\\ I_k(t) &= \varepsilon \int E_k(t)dt - \gamma \int I_k(t)dt\\\\ I_k(t) &= \varepsilon \int E_k(t)dt - R_k(t) \qquad \longrightarrow \varepsilon \int E_k(t) = I_k(t)+R_k(t) \end{align*}
Finally, the equation I cannot solve: \begin{align*} \int \frac{d}{dt} E_k(t) &= \int \lambda kS_k(t) \Theta(t) - \varepsilon E_k(t)\\\\ \int dE_k(t) &= \int \big( \lambda kS_k(t) \Theta(t) - \varepsilon E_k(t)\big) dt\\\\ E_k(t) &= \int \lambda kS_k(t) \Theta(t)dt - \varepsilon \int E_k(t)dt\\\\ E_k(t) + \varepsilon \int E_k(t)dt &= \lambda k \int S_k(t) \Theta(t)dt\\\\ E_k(t) + I_k(t) + R_k(t) &= \lambda k \int S_k(t) \Theta(t)dt\\\\ 1-S_k(t) &= \lambda k \int S_k(t) \Theta(t)dt \end{align*}
I do not know how to proceed from here. I was thinking if I could just have the integral $\int \Theta(t)dt$ then I could use the auxiliary function $\phi(t)$ but I do not know how.