I don't know the answer to this question, but the algebra of differential operators is almost commutative. So by looking at the principal symbols (with respect to the standard filtration by the degree of $D=\partial_x$), you can conclude that $\sigma(L_1)=\sigma(L_2)$. You can push it further a bit by considering other filtrations if, say, $L_i$ are polynomial coefficient operators, and by looking at the quasi-classical approximation given by the Poisson bracket.
NB: We assume that the coefficients of the differential operators belong to a commutative domain (polynomial, rational, algebraic, exp-algebraic, etc functions). The situation is quite different for various matrix-valued operators, where it is more reasonable to expect braid relations.
As an illustration, it is straightforward to show that in degree one, $L_1=L_2$. Indeed, it is well-known (and easy to prove) that any first order linear differential operator can be conjugated to $a\partial$ in a suitable differential extension (this requires adding exponentials of antiderivatives, as in solving $(\partial+p)f=0$). So without loss of generality, we may assume that $L_1=a\partial$ and $L_2=a\partial+b$, where $a,b$ are functions. Looking at the subprincipal terms, i.e. equating the coefficients of $\partial^2$ in the braid relation, we conclude that $a^2b=0$ and so $b=0$, i.e. $L_1=L_2$.