Completely ignoring all convergence issues, this really is just following your nose.
Plugging in the representation you give, you want to check
$$ \sum_{n = 1}^\infty (1 - e^{-t})^{n-1} Q_n^+(f_0) \overset{?}{=} f_0 + \int_0^t e^{s} Q^+(f_s,f_s) ~ds $$
Using the representation you give again to replace $f_s$, you rewrite this as checking (using that $Q^+$ is bilinear)
$$ \sum_{n = 1}^\infty (1-e^{-t})^{n-1}Q_n^+(f_0) \overset{?}{=} \\Q_1^+(f_0) + \sum_{n = 1}^\infty \sum_{m = 1}^\infty \int_0^t e^{-s} (1 - e^{-s})^{n-1} (1-e^{-s})^{m-1} Q^+(Q_n^+(f_0), Q_m^+(f_0)) ~ds $$
The integral can now be explicitly evaluated, as $\frac{d}{ds} (1-e^{-s})^k = k (1 - e^{-s})^{k-1} e^{-s}$, so we reduce to
$$ \sum_{n = 1}^\infty (1-e^{-t})^{n-1}Q_n^+(f_0) \overset{?}{=} Q_1^+(f_0) + \sum_{n,m = 1}^\infty \frac{1}{n+m-1} (1 - e^{-t})^{n+m-1} Q^+(Q_n^+(f_0), Q_m^+(f_0)) $$
Rewrite the sum in terms of $n$ and $k = n+m$, you get
$$ \sum_{n = 1}^\infty (1-e^{-t})^{n-1}Q_n^+(f_0) \overset{?}{=} Q_1^+(f_0) + \sum_{k = 2}^\infty (1 - e^{-t})^{k-1}\underbrace{\frac{1}{k-1} \sum_{n = 1}^{k-1} Q^+(Q_n^+(f_0), Q_{k-n}m^+(f_0))}_{= Q_k^+(f_0)} $$
and so equality follows.
Convergence can be checked for small $t$: if the bilinear form satisfies the bound $|Q^+(f,g)| \leq C|f||g| $, then by induction you can show that $|Q_n^+(f_0)| \leq C^{n-1} |f_0|^n$. For $|t|\ll 1$ you have that $|(1-e^{-t})| \ll 1$ and hence the representation converges as it is dominated by a geometric series. And so for $t$ sufficiently small all of the interchanges of integrals and sums, and reindexing of summation, are valid in the formal computation above.
This is more a long comment than an answer, since there's a peculiarity of definition \eqref{2} that I don't understand and makes me think the two definition aren't equivalent. To start, let's consider the part containing explicitly the time derivative in definition \eqref{2}: we have
$$
\begin{split}
\frac{\mathrm d}{\mathrm d t} \int_\Omega \phi_t ( x) \, \mathrm d \mu_t (x) &= \lim_{h\to 0} \frac{1}{h} \int_\Omega \big[\phi_{t+h} ( x)\, \mathrm d \mu_{t+h} (x) - \phi_t ( x) \, \mathrm d \mu_t (x)\big]\\
&= \lim_{h\to 0} \frac{1}{h} \int_\Omega \big[\phi_{t+h} ( x)\, \mathrm d \mu_{t+h} (x) -\phi_{t+h} ( x)\,\mathrm d\mu_t (x) \\
&\qquad\qquad\quad\quad +\phi_{t+h} ( x)\, \mathrm d\mu_t (x) - \phi_t ( x) \, \mathrm d \mu_t (x)\big]\\
&=\lim_{h\to 0} \frac{1}{h} \int_\Omega \phi_{t+h} ( x)\big[\mathrm d \mu_{t+h} (x) -\mathrm d\mu_t (x)\big] \\
&\qquad\qquad + \lim_{h\to 0} \frac{1}{h} \int_\Omega\Big[\phi_{t+h} ( x) - \phi_t ( x)\big] \mathrm d \mu_t (x)\\
&=\left.\frac{\mathrm d}{\mathrm d t} \int_\Omega \phi_{s} ( x)\,\mathrm d \mu_{t} (x) \right|_{s=t} + \left. \int_\Omega \frac{\mathrm d}{\mathrm d s}\phi_{s} ( x) \, \mathrm d \mu_t (x)\right|_{s=t}\\
&=\left.\frac{\mathrm d}{\mathrm d t} \int_\Omega \phi_{s} ( x)\,\mathrm d \mu_{t} (x) \right|_{s=t} + \int_\Omega \Big(\frac{\mathrm d}{\mathrm d t}\phi_{t} ( x) \Big) \, \mathrm d \mu_t (x)
\end{split}
$$
Loosely speaking, I got this identity by "freezing" the time dependence of the test function in the first term of the equation: this implies that if $\mu_t(x)$ is a solution of \eqref{a} according to definition \eqref{3}, then it is not a solution according to definition \eqref{2}, unless
we have
$$
\int_\Omega \Big(\frac{\mathrm d}{\mathrm d t}\phi_{t} ( x) \Big) \, \mathrm d \mu_t (x) \equiv 0 \quad\forall\phi \in C^\infty_c ((0, 1) \times \Omega)\label{4}\tag{vc}
$$
But this leads straight to the reason for my doubt: the vanishing condition \eqref{4} cannot be satisfied unless $\mu_t(x)\equiv0$ identically!
This can be seen by constructing a class of test functions $\{_{\tiny M,N}\phi_t(x)\}_{MN\in\Bbb N}$ by using any smooth partition of unity of the domain $\Omega$, say $\{\psi_k(x)\}_{k\in \Bbb N}$ and any smooth partition of unity of the closed interval $[0,1]$ say $\{\chi_n(t)\}_{n\in \Bbb N}$. Formally speaking,
$$
{_{\tiny M,N}\phi}_t(x) = t\sum_{n=1}^N\chi_n(t)\sum_{k=1}^M\psi_k(x) \quad M, N\in\Bbb N.
$$
Am I wrong?
Best Answer
Let's focus on $N=1$ only, the case $N>1$ is just a tensorization of the argument below.
The whole argument is actually unrelated to the specific (aggregation-diffusion) PDE or gradient flow: As soon as you have a curve of probability measures $\mu_t$ satisfying the continuity equation $$ \partial_t\mu_t+\nabla\cdot(v_t\mu_t)=0 $$ for some (smooth enough) vector-field $v=v_t(x)$, then in fact you have the representation $$ \mu_t=\Phi_t\#\mu_0, $$ where $\Phi_t(x)$ is the flow-map associated to the ODE $\begin{cases}\dot X_t=v_t(X_t)\\ X_0=x\end{cases}$. This is easy to check formally, and for a precise reference I can suggest reading chapter 8 in [1] (e.g. lemma 8.1.6 or prop. 8.1.8).
Now, for your specific context, it is just a matter of rewriting the PDE into a continuity equation. Factoring out and explointing the classical identity $\Delta\mu=\nabla\cdot\nabla\mu=\nabla\cdot(\frac{\nabla\mu}{\mu}\mu)=\nabla\cdot((\nabla\log \mu) \mu)$ you see that $\mu$ solves the continuity equation $$ \partial_t\mu_t=\nabla\cdot(\nabla\mu_t-(b\ast\mu_t)\mu_t) \quad\Leftrightarrow \quad \partial_t\mu_t+\nabla\Big((\underbrace{b\ast\mu_t-\nabla\log\mu_t}_{v_t})\mu_t\Big)=0. $$ Then the flow-map $\Phi_t$ can be simply represented from the implicit Duhamel formula as $$ \begin{cases}\dot X_t=v_t(X_t)\\ X_0=x\end{cases} \quad\Leftrightarrow\quad X_t=x+\int_0^t v_s(X_s)ds. $$ Plugging in the explicit formula for the velocity field $v_t=b\ast\mu_t-\nabla\log\mu_t$ gives exactly your expression (2).
[1] Ambrosio, L., Gigli, N., & Savaré, G. (2008). Gradient flows: in metric spaces and in the space of probability measures. Springer Science & Business Media.