The differential equation $\dot{\mathbf{x}}(t)=M(t)\mathbf{x}(t)$, conditions for a specific type of matrix $M$ to commute at different times I have some $N$x$N$ matrices $M_N(t)$ of the general form
$$\tag{1}
M_4(t)=\left[\begin{matrix}-f_0 & f_1 & 0 & 0  \\
 0 & -f_1 & f_2 & 0  \\
0 & 0 & -f_2 & f_3  \\
0 & 0 & 0 & -f_3 
\end{matrix}\right] \qquad M_5(t)=\left[\begin{matrix}-f_0 & f_1 & 0 & 0  & 0\\
 0 & -f_1 & f_2 & 0 & 0 \\
0 & 0 & -f_2 & f_3 & 0 \\
0 & 0 & 0 & -f_3 & f_4 \\
0 & 0 & 0 &  0 &-f_4
\end{matrix}\right]
$$
where each $f_n=f_n(t)$. From the physics of the problem, none of the $f_n$ may vanish for all $t$, and they are all non-negative. These matrices arise in studying the differential equation
$$\tag{2}
\dot{\mathbf{x}}(t)=M(t)\mathbf{x}(t)
$$
where $\mathbf{x}(t)$ is an $N$-tuple $(x_0(t), \ x_1(t), \ \dots, x_{N-1}(t))^\top$. I would like to study the conditions on $f_n$ under which the matrices $M$ commute at unequal times, that is
$$\tag{3}
\left[ \int\limits_0^t dt' M(t'), \ M(t)  \right]\stackrel{?}{=}0
$$
where $[ \cdot , \cdot]$ is the commutator. This is important to me because under these conditions the solution to (2) is simply $\mathbf{x}(t)=\exp\left( \int\limits_0^t dt' M(t')\right)\mathbf{x}(0)$ while without condition (3) the solution is a time ordered exponential.
By playing around with small $N$ cases in a CAS, I think that condition (3) may only be fulfilled only if
$$\tag{4}
M(t)=g(t) M_0
$$
where $g(t)$ is a function and $M_0$ is a constant matrix. Obviously (4) is a sufficient condition for (3), but I am uncertain about its necessity. To demonstrate this is a nontrivial claim, note that if we allow any of the $f_n$ to vanish for all $t$, then $M$ may be partly block diagonalized and (4) is not a necessary condition.
Question: Is (4) a necessary condition for (3) to hold with $f_n(t)\neq 0$, or is there another condition?
Related question: given the specific form of $M$, can the general time ordered solution $\mathbf{x}(t)=\mathcal{T}\exp\left( \int\limits_0^t dt' M(t')\right)\mathbf{x}(0)$ be simplified at all? Here $\mathcal{T}$ denotes the time ordering symbol.
Background (not necessary for the question)
Equation (2) arises from a master equation governing the number of electrons trapped in a quantum dot. The functions $f_n(t)$ are complicated, but the controlling factor of each is of the form $f_n(t) \approx \exp(e^{\alpha_n t})$
 A: Your sufficient condition appears to also be necessary.

Let $h_n(t)=\int_0^t f_n(\tau)d\tau$, and $H_n(t)=\int_0^t M_n(\tau)d\tau$.
Play around with small $n$ cases, e.g:
$$[M_2,H_2] = \begin{bmatrix} 0 & h_0f_1-h_1f_0\\ 0 & 0\end{bmatrix}$$
Look at $(j+2)^\text{th}$ row and $(j+1)^\text{th}$ column of $[M_n,H_n]$ to find:
$$ h_jf_{j+1}=h_{j+1}f_j \quad \forall_{j=0,\cdots,n-2}$$
Note that this implies $f_{j+1}(t)\propto f_j(t)$ since:
$$\require{cancel}\begin{align}
&h_jf_{j+1}=h_{j+1}f_j \\
\Leftrightarrow &f_{j+1}\int_0^t f_j = f_{j}\int_0^t f_{j+1} \tag{1} \\
\stackrel{d/dt}{\Rightarrow} &\bcancel{f_{j+1}f_{j}} + f'_{j+1}\int_0^t f_{j} = \bcancel{f_{j+1}f_{j}} + f'_{j}\int_0^t f_{j+1} \tag{2}
\end{align} $$
$$\begin{align}
(1) \Rightarrow f_{j+1}(t) &= \alpha_{j(j+1)}(t)f_j(t)\\
\stackrel{d/dt}{\Rightarrow} f'_{j+1}(t) &= \alpha'_{j(j+1)}(t)f_j(t)+\alpha_{j(j+1)}(t)f_j'(t) \tag{4}\\
(2) \Rightarrow f'_{j+1}(t) &= \alpha_{j(j+1)}(t)f'_j(t)\tag{5}
\end{align} $$
Together with $f_j(t)\neq 0$, $(4)$ and $(5)$ imply $\alpha_{j(j+1)}'(t)=0$.
Summarizing, we found that $\begin{cases}[M_n,H_n]=0\\ f_j(t)\neq 0 \quad \forall_{j=0,\cdots,n-1}\end{cases}$ if and only if $f_{j+1}(t)= \alpha_{j(j+1)} f_j(t)\ \forall_{j=0,\cdots,n-2}$, which can be recast as
$$f_j(t) = \beta_j f_0(t) \quad \forall_{j=1,\cdots,n-1}$$
