Matrix exponential - general form of the solution Motivation
This question is related to continuous-time Markov chains and models of DNA evolution. The question is asked in complete generality, however.
Background
Let $\boldsymbol Q$ be a $4 \times 4$  such that 
$$ \boldsymbol Q = \begin{bmatrix}
-( ap_2 + bp_3 + cp_4) & ap_2 & bp_3 & cp_4\\
ap_1 & -(ap_1 + dp_3 + ep_4 )  & dp_3& ep_4  \\
bp_1 & dp_2 & -( bp_1 + dp_2 + p_4 )& p_4  \\
cp_1 & ep_2 & p_3 & -(cp_1 + ep_2 + p_3)  \\
\end{bmatrix} $$
with all off-diagonal entries are positive, with rates $a \ldots e > 0$ and equilibrium probabilities $ 0 < p_i < 1, \, i = 1, 2, 3, 4$ and $\sum_i p_i = 1$.
For $t > 0$, we can then construct the stochastic matrix $\boldsymbol P(t) := \exp(t\boldsymbol Q)$. 
Example: as a special case, consider $p_i = 1/4$ for all $i$ and $a = b = \ldots = e = \lambda$. Then $\boldsymbol P(t)$ is of the form
$$ \boldsymbol P(t) = \begin{bmatrix}
p_0(t) & p_1(t) & p_1(t) & p_1(t)\\
p_1(t) &  p_0(t) & p_1(t)&  p_1(t)  \\
p_1(t) & p_1(t) & p_0(t)& p_1(t)  \\
p_1(t) & p_1(t) &p_1(t) &  p_0(t) \\
\end{bmatrix} $$
whence $p_0(t) = \frac{1}{4} + \frac{3}{4}\exp(-4\lambda t)$ and $p_1(t)=   \frac{1}{4} - \frac{1}{4}\exp(-4\lambda t)$.
Question
I am aware it is not possible to obtain $\boldsymbol P(t)$ in closed-form for the general formulation. 
For computation, diagonalisation is the way to go. 

My question is: is it possible to at least know the general functional form the entries in $\boldsymbol P(t)$ will take?

For instance, in the example above the general form of the solution(s) is $ p(t) = \alpha + \beta \exp(-\gamma t)$.
For more general models the pattern seems to sort of carry on, with solutions looking like $p(t) = \sum_{j = 1}^K w_j\exp(-a_j t)$. Are there any tools/(simple) facts I can use to prove this, should it be true?
 A: Indeed, the entries in the exponential $P(t)=e^{tQ}$ will be of the form $\sum w_j\exp(a_jt)$ (i.e. linear combinations of exponentials), so long as the matrix $Q$ is diagonalizable. In this case, there exists an invertible matrix $B\in M_{4\times 4}(\mathbb R)$ such that
$$B^{-1}QB=D=\mathrm{Diag}\{\lambda_1,\lambda_2,\lambda_3,\lambda_4\}$$
and therefore 
$$(tQ)=B(tD)B^{-1}=B\cdot \mathrm{Diag}\{t\lambda_1,t\lambda_2,t\lambda_3,t\lambda_4\}\cdot B^{-1}.$$
As you may know, the exponential can then be calculated by the power series:
$$e^{tQ}=\sum_{k=0}^\infty \dfrac{t^kQ^k}{k!} =B\left( \sum_{k=0}^\infty \dfrac{t^kD^k}{k!} \right)B^{-1}=B\cdot \mathrm{Diag}\{e^{t\lambda_1},e^{t\lambda_2},e^{t\lambda_3},e^{t\lambda_4}\}\cdot B^{-1}.$$
Therefore, if $Q$ is diagonalizable, the entries of $e^{tQ}=(a_{ij}(t))
$ are linear combinations of the form
$$a_{ij}(t)=w_{ij,1}\,e^{t\lambda_1}+w_{ij,2}\,e^{t\lambda_2}+w_{ij,3}\,e^{t\lambda_3}+w_{ij,4}\,e^{t\lambda_4},$$
for some $w_{ij,k}$ ($k=1,2,3,4$) which are determined by $B$ and $B^{-1}$. If one of the eigenvalues $\lambda_i$ is $0$ you might get a constant term like in your example.
However, real matrices are not always diagonalizable -- in the language of topology, they're not a dense set in $M_{n\times n}(\mathbb R)$, which means that in the set of matrices $Q$ given by your definition (with $9$ parameters $p_1,...,p_4$ and $a,...,e$) you're very likely to stumble upon a non-diagonalizable matrix, even if it's an stochastic matrix.
Luckily, in the case where $Q$ is not diagonalizable we can take the real Jordan form, which always exists. That is, we can always find an invertible matrix $B$ such that $B^{-1}QB=J$, where $J$ is a real, block diagonal matrix of the form
$$J=\begin{pmatrix}J_1\end{pmatrix}, \begin{pmatrix}J_1 & 0 \\ 0 & J_2\end{pmatrix}, \begin{pmatrix}J_1 & 0 & 0 \\ 0 & J_2 & 0 \\ 0 & 0 & J_3\end{pmatrix} \text{or} \begin{pmatrix}J_1 & 0 & 0 & 0 \\ 0 & J_2 & 0 & 0 \\ 0 & 0 & J_3 & 0 \\ 0 & 0 & 0 & J_4\end{pmatrix}$$
where each $J_i$ is a square matrix of size $1\times 1, ..., 4\times 4$. (If there are four $1\times 1$ blocks then $J$ is a diagonal matrix and $Q$ is diagonalizable.)
Furthermore, the exponential of a Jordan form is very well known, so calculating $e^{tQ}=Be^{tJ}B^{-1}$ is relatively easy. Since $J$, and therefore $tJ$ is block diagonal, the exponential $e^{tJ}$ is also block diagonal, where the blocks are the corresponding exponentials of the blocks in $tJ$.
Since listing all the blocks and the respective exponentials that could possibly appear would lengthen this already long answer considerably (and is a topic already covered in countless sources), I'll leave you with the following result from Lawrence Perko's Differential equations and dynamical systems (2006, p.42):

Corollary. Each coordinate in the solution $x(t)$ of the initial value problem $x'=Ax$ is a linear combination of functions of the form $$t^ke^{at}\cos bt ~~~\text{or}~~~ t^ke^{at}\sin bt$$ where $\lambda=a+ib$ is an eigenvalue of the matrix $A$ and $0\leq k \leq n-1$.

Since here $x(t)=e^{At}x_0$ for some initial condition $x_0\in \mathbb R^n$, this applies directly to your problem: we can conclude that the entries of the exponential $e^{Qt}$, where $Q$ is defined as in your question, are linear combinations of terms of the form
$$e^{a_kt}\cos b_kt, e^{a_kt}\sin b_kt, te^{a_kt}\cos b_kt, ..., t^3e^{a_kt}\sin b_kt$$
where $\lambda_k= a_k+ib_k$ is one of the eigenvalues of $Q$ ($k=1,2,3,4$ with some $\lambda_k$ possibly equal, i.e. of multiplicity $\geq 2$). Moreover, if $Q$ is an stochastic matrix, then each of its eigenvalues has a norm lesser than or equal to $1$, and $1$ is always an eigenvalue of $Q$. I believe this is as far as you can get without imposing excessive restrictions on the parameters $p_1,...,p_4$ and $a,...,e$.
For further reading on the Jordan form and matrix exponentials, section 1.8 of Perko's book is a good source, but it's also covered in some linear algebra and differential equations textbooks (particularly those that treat linear systems like Perko does).
