Epidemic threshold on a directed network

I am trying to understand the equations used in a paper

Mainly I'm trying to understand how the epidemic thershold was calculated using the equation: $$I^{t+\Delta t}_a = -\mu\Delta tI^{t}_a + I^{t}_a + \lambda(N^{t}_a - I^{t}_a)a\Delta t\int d{a}'\frac{I^{t}_a}{N}+\lambda(N^{t}_a - I^{t}_a)\int d{a}'\frac{I^{t}_a{a}'\Delta t}{N}$$

From the supplementary information (https://media.nature.com/original/nature-assets/srep/2012/120625/srep00469/extref/srep00469-s1.pdf (eq.17)), I get that $$I^{t}_{a} - \text{number of infected individuals in the class a, at time t}$$ $$N^{t}_a - \text{total number of individuals with activity a, at time t}$$ $$\lambda - \text{transition probability, } \mu - \text{recovery probability}$$

I understand that first two terms $-\mu\Delta tI^{t}_{a}$ and $I^{t}_{a}$ represent the number of recovered individuals, and number of infected individuals of class $a$ at time $t$ respectively. Rest of the equation is a mystery however.

In the supplementary information eq.18, we get: $$\int daI^{t+\Delta t}_a = I^{t + \Delta t} = I^{t} - \mu\Delta tI^{t} + \lambda \langle a \rangle I^{t}\Delta t + \lambda \theta^{t}\Delta t$$

where $\theta^{t} = \int d{a}'I^{t}_{{a}'}{a}'$, I don't understand why we get rid of N in the $4^{th}$ term and how we obtain a moment of a distribution from the $2^{nd}$ term to get this. We then get:

$$\theta^{t+\Delta t} = \theta^{t} - \mu\theta^{t}\Delta t + \lambda \langle a^{2} \rangle I^{t}\Delta t + \lambda \langle a \rangle \theta^{t}\Delta t$$

Which again makes no sense to me, I feel like there is a lot of steps missing here. From here onward it's just dark magic for me, no clue what's going on.

Then this is written as differential forms: $$\partial_tI = -\mu I + \lambda \langle a \rangle I + \lambda\theta$$ $$\partial_t\theta = -\mu\theta + \lambda \langle a^{2} \rangle I + \lambda \langle a \rangle \theta$$ From that somehow the Jacobian is calcuated as $$J = \begin{pmatrix} -\mu + \lambda \langle a \rangle & \lambda \\ \lambda \langle a^{2} \rangle & -\mu + \lambda \langle a \rangle \end{pmatrix}$$

Eigenvalues: $$\Lambda_{(1,2)} = \lambda \langle a \rangle - \mu \pm \lambda \sqrt{\langle a^{2} \rangle}$$

Threshold: $$\frac{\lambda}{\mu} > \frac{1}{\langle a \rangle + \sqrt{\langle a^{2} \rangle}} + \mathcal{O}(\frac{1}{N})$$

I would appreciate if someone could explain this to me, I have been trying for few weeks now, and it just makes less sense more I look at it.

Now that's the threshold on an undirected network, if I were to calculate the threshold on a directed network how would I do that? Would it be something like this:

$$d_{t}I = -\mu Ia + \lambda(N_{a}-I_{a})a\sum_{{a}'}I_{{a'}}\frac{m}{N} + (N_{a} - I{a})\lambda \sum_{{a'}}{a}'I_{{a}'}\frac{m}{N}$$