Setup
We have a system of $N$ diffusion processes described by
$$X_t^{i}=X_0^{i}+\int_0^t \mu(s,X_s)ds+\int_0^t\sigma(s,X_s)dW^{i}_s,$$
with the Brownian motions $W^{i}$ as well as the drift $\mu$ and the volatility $\sigma$ with sufficient regularity (e.g. Lipschitz condition and sublinear growth) for strong uniqueness and existence of a solution.
On this basis, define the process $$\Lambda_t^{i}=|X_t^{i}|+1/N\sum_{j=1}^N|X_t^{j}|,$$
and assume it to be Sub-Gaussian distributed (i.e. we can find a Gaussian random variable as an upper bound).
Problem
In this setup I'd like to get the estimation
\begin{equation}\label{thing}
\mathbb{E}[(1+\Lambda_i^q)\mid X_0^{i}=x]\leq C(1+x^q)\quad \text{where }C>0\text{, }q>1\text{ and }x\in\mathbb{R}.
\end{equation}
Idea
The basic idea for proving this was using Gronwall's inequality in the monotonic case. So if I could estimate
$$\Lambda_t^{i}\leq\alpha(t)+\int_0^t\beta(s)\Lambda_s^{i}ds,$$
for $\alpha$ monotonoically increasing, that would give me
$$\Lambda_t^{i}\leq \alpha(t)e^{\int_0^t\beta(s)ds}.$$
My approach would be to choose $\alpha(t)=X_0^{i}$ or $\alpha(t)=X_0^{i}+1/N\sum_{j=1}^N X_0^{j}$ for some $\beta$ I'm not able to find. But it would yield
\begin{align}
\mathbb{E}[1+(\Lambda_t^{i})^q)|X_0^{i}=x]\leq (1+Cx^q)
\end{align}
for some constant, since sub-gaussianity of $\Lambda_t^{i}$ implies sub-gaussianity of the $X_t^j$ represented by this constant $C$.
Question
Do you have any idea how to use Gronwall's inequality in this case, how to choose $\beta$ or $\alpha$ to get to
\begin{equation}
\mathbb{E}[(1+\Lambda_i^q)\mid X_0^{i}=x]\leq (1+Cx^q)\quad \text{where }C>0\text{, }q>1\text{ and }x\in\mathbb{R}.
\end{equation}
Edit: For $\epsilon>0$ sufficiently small, one can use the sub-gaussianity to follow $\mathbb{E}[e^{\epsilon \sup_{0\leq t\leq T} (\Lambda_t^{i})^2}]\leq C_{T,\epsilon}$, on a time horizon $[0,T]$ and $C_{T,\epsilon}>0$ depending on $T$ and $\epsilon$. Maybe this can be of help.