See EDIT #1 for numerical Solutions.
As a first step I would try to simplify the original system of non-linear ODEs, solve the simplified version, and look for the asymptotic result.
The equations are
eq0 = {
l'[t] == sl - ml*(M[t]*l[t])/(1 + l[t]) - dl*l[t],
M'[t] ==
sM (a1*l[t]/(B + l[t]) + a2*(M[t]*l[t])/(1 + l[t])) -
vM*(M[t]*h[t])/(1 + h[t]),
h'[t] == sh - vh*(M[t]*h[t])/(1 + h[t]) - dh*h[t]
};
eq0 can't be solved by DSolve[]:
sol = DSolve[eq0, {l[t], M[t], h[t]}, t] // Simplify
(*
Out[19]= DSolve[{l[t] (dl + (ml M[t])/(1 + l[t])) + Derivative[1][l][t] ==
sl, (vM h[t] M[t])/(1 + h[t]) + Derivative[1][M][t] ==
l[t] ((a1 sM)/(B + l[t]) + (a2 sM M[t])/(1 + l[t])),
h[t] (dh + (vh M[t])/(1 + h[t])) + Derivative[1][h][t] == sh}, {l[t], M[t], h[t]}, t]
*)
First attempt: linearization
We set all non-linear terms to zero.
This gives
eqlin = {
l'[t] == sl - dl*l[t],
M'[t] == sM (a1*l[t]/B),
h'[t] == sh - dh*h[t]
};
The solution is
sollin = DSolve[eqlin, {l[t], M[t], h[t]}, t] // Simplify
(*
Out[15]= {{
l[t] -> sl/dl + E^(-dl t) C[1],
M[t] -> -((a1 sM (sl - dl sl t + dl (-1 + E^(-dl t)) C[1]) - B dl^2 C[2])/(B dl^2)),
h[t] -> sh/dh + E^(-dh t) C[3]}}
*)
We see that the solution for t->[Infinity] is, provided dl > 0:
sollin /. {C[3] -> 0, C[2] -> 0, C[1] -> 0}
(*
Out[28]= {{
l[t] -> sl/dl,
M[t] -> -((a1 sM (sl - dl sl t))/(B dl^2)),
h[t] -> sh/dh}}
*)
Which means that l[t] and h[t] remain finite while M[t] increases indefinitely with t.
Second attempt: large functions l[t] and h[t]
In the limit l[t]->[Infinity], h[t]->[Infinity] the system becomes
eqinf = {
l'[t] == sl - ml*(M[t]) - dl*l[t],
M'[t] == sM (a1*1) + a2*M[t] - vM*M[t],
h'[t] == sh - vh*(M[t]) - dh*h[t]
};
It is linear and the solution is
solinf = DSolve[eqinf, {l[t], M[t], h[t]}, t] // Expand
(*
Out[30]= {{
h[t] -> sh/(a2 + dh - vM) + (a2 sh)/(dh (a2 + dh - vM)) + (a1 sM vh)/(
dh (a2 + dh - vM)) + (
a1 E^(-t (a2 + dh - vM)) sM vh)/((a2 - vM) (a2 + dh - vM)) - (sh vM)/(
dh (a2 + dh - vM)) - (a1 sM vh)/((a2 + dh - vM) (-a2 + vM)) + (
a1 E^(-dh t - t (a2 - vM)) sM vh)/((a2 + dh - vM) (-a2 + vM)) +
E^(-dh t) C[1] + (E^(-dh t) vh C[3])/(a2 + dh - vM) - (
E^(t (a2 - vM)) vh C[3])/(a2 + dh - vM),
l[t] -> sl/(a2 + dl - vM) + (a2 sl)/(dl (a2 + dl - vM)) + (a1 ml sM)/(
dl (a2 + dl - vM)) + (
a1 E^(-t (a2 + dl - vM)) ml sM)/((a2 - vM) (a2 + dl - vM)) - (sl vM)/(
dl (a2 + dl - vM)) - (a1 ml sM)/((a2 + dl - vM) (-a2 + vM)) + (
a1 E^(-dl t - t (a2 - vM)) ml sM)/((a2 + dl - vM) (-a2 + vM)) +
E^(-dl t) C[2] + (E^(-dl t) ml C[3])/(a2 + dl - vM) - (
E^(t (a2 - vM)) ml C[3])/(a2 + dl - vM),
M[t] -> (a1 sM)/(-a2 + vM) + E^(t (a2 - vM)) C[3]}}
*)
The behaviour for large t should be considered carefully. We assume that all exponentials should go to zero for t->[Infinity].
This requires the following relation to hold.
a2 < vM < a2 + dh
Then the asymptotic solution is after some simplifications
solasy = {
l[t] -> sl/dl - (a1 ml sM)/(dl (vM - a2)),
M[t] -> a1 sM,
h[t] -> sh/dh - (a1 vh sM)/(dh (vM - a2)),
}
We see that the solution should be good for large enough sh and sl. The first terms of l[t] and h[t] coincide with those of the lienearized equations.
So at least we have found the asymptotic solution in two simplified cases. In the second case these solutions are even stationary. This might be a start.
EDIT #1
In order to explore the types of the full time dependent solutions of the ODEs, we present a code which can be used to experiment. We can also cautiously draw some first conclusions which includes a qualitative comparision of the solution types with the results of the simplifying procedures given earlier and kind of "backward justification" of the latter.
The code solves the time dependent problem numerically.
To be able to do this, we need values for the parameters and for the initial conditions. All that was gives in the OP is that all quantities are not negative. Hence we shall chose all these unknown quantities from a certain simple random source.
Preparation
The system of equations
eq0 = {
l'[t] == sl - ml*(M[t]*l[t])/(1 + l[t]) - dl*l[t],
M'[t] ==
sM (a1*l[t]/(B + l[t]) + a2*(M[t]*l[t])/(1 + l[t])) -
vM*(M[t]*h[t])/(1 + h[t]),
h'[t] == sh - vh*(M[t]*h[t])/(1 + h[t]) - dh*h[t]
};
The random source (giving rational numbers >0)
r := Rationalize[0.1 + 3*Random[], 1/10]
The source for the 11 parameters
par = {sl, ml, dl, sM, a1, B, a2, vM, sh, vh, dh};
is
r11 := Array[r &, 11];
The source for the 3 initial values
iv = {l0, M0, h0};
is
r2 := Array[r &, 3]
Running code
To be run "by hand".
Instruction to experiment:
First run the three blocks in the given order.
Then you can choose to change only the parameters leaving the inital values unchanged or vice versa, or choose to change both.
(* choose inital values *)
ri = r2;
repi = Thread[iv -> ri];
repis = ToString[repi, InputForm];
eq0ri = l[0] == l0 && M[0] == M0 && h[0] == h0 /.
repi(* equations with parameters and initials values set *);
(* choose parameters *)
rp = r11;
rep = Thread[par -> rp]; (* eins *)
reps = ToString[Take[rep, {1, 5}], InputForm] <> "\n" <>
ToString[Take[rep, {6, 11}], InputForm];
eq0r = eq0 /. rep (* equations with parameters set *);
(* solve ODEs in the time range from 0 to tmax, and display solution *)
tmaxc = 10; (* set maximum calculation time *)
tmaxd = 10; (* set maximum display time *)
sol = NDSolve[eq0r && eq0ri, {l[t], M[t], h[t]}, {t, 0, tmaxc}][[1]];
s = {l[t], M[t], h[t]} /. sol;
Plot[s, {t, 0, tmaxd}, AxesLabel -> {"t", "f[t]"},
PlotLabel ->
"Solution of ODEs boschbird\nparameters = " <> reps <>
"\ninitial values = " <> repis <>
"\ncurves : blue = l[t], yellow = M[t],green = h[t]"]
(* 150925_Plot_ODE_boschbird_1.jpg *)

Some conclusions
a)
We find solutions of two types
Type 1: l[t] and h[t] have finite asymptotic values, while M increases indefinitely linearly with time
Type 2: all three quantities l[t], M[t], h[t] have finite asymptotic values
In the range of random varibales chosen here, the majority of the solutons is of type 1.
b)
Type 1 solutions resemble qualitatively the solutons of simplification attempt 1
Type 2 solutions resemble qualitatively the solutons of simplification attempt 2
Closer inspection could reveal more than resemblance.
Final remarks
The code can be improved, of course, in many respects. For instance
The procedure of an overall random source can be replaced by sources adapted more closely to each parameter in question.
Once the numerical value of a parameter or an initial value is know, the random choice of it can be replaced.