Finding equilibria to a system of 3 ODEs with 11 unknown parameters. Looking for the intersection of 3 surfaces

Question

I have the following system of 3 ODEs of 3 variables $(\ell,M,h)$ and 11 parameters $(\sigma_{\ell},\mu_{\ell},d_{\ell},\sigma_M,\alpha_1,\beta,\alpha_2,\nu_M,\sigma_h,\nu_h,d_h)$:

$\frac{d\ell}{dt}=\sigma_{\ell} - \mu_{\ell}\frac{M \ell}{1+\ell} - d_{\ell} \ell$,
$\frac{dM}{dt}=\sigma_M\left(\alpha_1\frac{\ell}{\beta+\ell} +\alpha_2\frac{M \ell}{1+\ell} \right) - \nu_M \frac{Mh}{1+h} $,
$\frac{dh}{dt}=\sigma_h-\nu_h \frac{Mh}{1+h} - d_h h$.

Naturally I tried to simultaneously solve the system of nullclines as follows:

Reduce[
  sl - ml*(M*l)/(1 + l) - dl*l == 0 && 
   sM (a1*l/(B + l) + a2*(M*l)/(1 + l)) - vM*(M*h)/(1 + h) == 0 && 
   sh - vh*(M*h)/(1 + h) - dh*h == 0 && l >= 0 && M >= 0 && h >= 0,
  {l, M, h},
  Reals
 ]

But my laptop and the university computers ran indefinitely and in both cases I aborted the evaluation after about 3 hours.

I am really stuck and have no Idea how I can range these parameter values until a particular combination of values presents itself representing an equilibria solution to my ODE system. Any advice would help.

P.S. I have tried nondimensionalising the system and I was only able to reduce the parameter dimension by 1. And note that the parameters and variables must be real and positive (or zero).

Answer 1

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.