Excitation loss can be modelled using a Lindblad equation (with a Hermitian Hamiltonian) which is guaranteed to be trace-preserving. However, one must work in a slightly larger Hilbert space. In particular, we need to introduce one more state, the vacuum state $\lvert \mathrm{vac} \rangle$, representing the situation where there is no excitation in the network. Excitation loss is then modelled by an incoherent process which transfers population to the vacuum state.
The non-Hermitian Hamiltonian $H_\mathrm{loss}$ should be replaced by a Lindblad generator
$$ \mathcal{L}_\mathrm{loss} (\rho) = \sum_m \Gamma_m \left ( L_m \rho L_m^\dagger - \frac{1}{2} \lbrace L^\dagger_m L_m,\rho \rbrace\right),$$
where $\rho$ is the density matrix and $L_m = \lvert \mathrm{vac}\rangle\langle m \rvert$. Each Lindblad operator $L_m$ represents incoherent loss at site $m$.
How do we prove that these formulations are equivalent? Let's first see how the non-Hermitian formulation works. The Schroedinger equation for the state ket $\lvert \psi\rangle$ is ($\hbar =1$)
$$ \frac{\mathrm{d}}{\mathrm{d}t}\lvert \psi\rangle = -\mathrm{i}(H_0 + H_\mathrm{loss}) \lvert \psi\rangle = \left(-\mathrm{i}H_0 - \sum_m \frac{\Gamma_m}{2} \Pi_m \right) \lvert \psi\rangle, $$
where $\Pi_m = \lvert m \rangle \langle m \rvert$. The evolution of the state bra $\langle \psi\rvert$ is given by the Hermitian conjugate of the above equation. From these one derives the von-Neumann equation for the density matrix
$$ \frac{\mathrm{d}\rho}{\mathrm{d}t} = -\mathrm{i}[H_0,\rho] - \sum_m\frac{\Gamma_m}{2}\{\Pi_m,\rho\}.$$
It is easier to work in the Heisenberg picture, where the equation of motion for an operator $A$ can be found using $$ \frac{\mathrm{d}\langle A\rangle}{\mathrm{d}t} = \mathrm{Tr} \left[ A \frac{\mathrm{d}\rho}{\mathrm{d}t} \right ] = \mathrm{Tr} \left[\frac{\mathrm{d}A}{\mathrm{d}t} \rho \right ]. $$
From this we deduce
$$ \frac{\mathrm{d}A}{\mathrm{d}t} = \mathrm{i}[H_0,A] - \sum_m\frac{\Gamma_m}{2}\{\Pi_m,A\}.$$
To see that the anti-commutator implies excitation loss, let's ignore $H_0$ for the moment (i.e. set $H_0 = 0$). Then the population of site $m$ is given by the solution of
$$ \frac{\mathrm{d}\langle \Pi_m\rangle}{\mathrm{d}t} = -\Gamma_m \langle \Pi_m\rangle,$$
or in other words $\langle \Pi_m\rangle(t) = \mathrm{e}^{-\Gamma_m t}\langle \Pi_m\rangle(0)$. The population of each site decays exponentially at a rate $\Gamma_m$. This obviously leads to a decay of the total probability (i.e. the trace of the density matrix) in time:
$$ \mathrm{Tr} [\rho(t)] = \sum_m \langle \Pi_m\rangle(t) = \sum_m \mathrm{e}^{-\Gamma_m t}\langle \Pi_m\rangle(0).$$
Now examine the Lindblad equation
$$ \frac{\mathrm{d}\rho}{\mathrm{d}t} = -\mathrm{i}[H_0,\rho] + \mathcal{L}_\mathrm{loss}(\rho).$$
In the Heisenberg picture, one finds
$$ \frac{\mathrm{d}A}{\mathrm{d}t} = \mathrm{i}[H_0,A] + \mathcal{L}_\mathrm{loss}^\dagger(A),$$
where the adjoint Liouvillian is
$$ \mathcal{L}^\dagger_\mathrm{loss} (A) = \sum_m \Gamma_m \left ( L_m^\dagger A L_m - \frac{1}{2} \lbrace L^\dagger_m L_m,A \rbrace\right).$$
This generates the same dynamics as the non-Hermitian formulation in the single-excitation subspace. (The single-excitation subspace is spanned by the set of states with one excitation, i.e. the states {$\lvert m \rangle\}$.) To prove this, we show that the Heisenberg equation for any population or coherence in the single excitation subspace (i.e. any operator $\lvert m\rangle \langle n\rvert$) is the same in either formulation. This follows because $L_k^\dagger \lvert m\rangle \langle n\rvert L_k = 0$ for all $k, m, n$. Any operator $A_1$ acting in the single-excitation subspace can be represented as a linear combination of the operators $\lvert m\rangle \langle n\rvert$, and therefore
$$ \frac{\mathrm{d}A_1}{\mathrm{d}t} = \mathrm{i}[H_0,A_1] - \sum_m\frac{\Gamma_m}{2}\{\Pi_m,A_1\}.$$
This proves that the two formulations give identical results for observables in the single-excitation subspace.
However, the equation of motion for the population of the vacuum state is
$$\frac{\mathrm{d}}{\mathrm{d}t} \langle \lvert \mathrm{vac} \rangle\langle \mathrm{vac}\rvert \rangle = \sum_m \Gamma_m \langle \Pi_m\rangle, $$
which shows that the single-excitation populations $\langle \Pi_m\rangle$ act as an (incoherent) source term for the vacuum population. This ensures preservation of the trace.
