1
$\begingroup$

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).

|improve this question
$\endgroup$
  • 1
    $\begingroup$ Out of curiosity, what does the system represent? $\endgroup$ – Chris K Sep 23 '15 at 19:02
  • 2
    $\begingroup$ This system represents the development of an Atherosclerotic plaque within the intima of an artery $\endgroup$ – boschbird Sep 24 '15 at 4:37
  • $\begingroup$ @ boschbird: very interesting application (and of practical importance for many people). Even without knowing this, I found the non-linear mathematical / Mathematica problem interesting enough to invest more time. You can find the result in EDIT #1 to my answer. $\endgroup$ – Dr. Wolfgang Hintze Sep 25 '15 at 9:09
  • $\begingroup$ @Dr.WolfgangHintze Thank you very much. I really appreciate this help. If you are still interested in the mathematical analysis of Atherosclerosis I've attached a link to a preprint of an article that my supervisor and another PhD collaborated Chalmers_etal_2015 $\endgroup$ – boschbird Nov 17 '15 at 0:19
  • $\begingroup$ @boschbird: Thank you for providing the link to the interesting paper. At a first glance I found the bifurcation most interesting. This, of course needs non-linear dynamics. Just in passing, I remember the Hagen-Poisseuille law from the time of my physics study. It shows the strong sensitivity of the flux on the radius (r^4). $\endgroup$ – Dr. Wolfgang Hintze Nov 17 '15 at 17:14
4
$\begingroup$

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 *)

enter image description here

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.

|improve this answer
$\endgroup$
  • $\begingroup$ I think that both of these approaches have serious flaws: 1) Linearization: what justification is there to simply 'drop' the non-linear terms. All linearizations of the system should involve Taylor expansions about an equilibria (i.e. we obtain the Jacobian matrix of the system by dropping higher order terms). Also, your application of your method of 'linearization' is incorrect; the 2nd equation should be "M'[t] == sM*a1" and in this way we are left with 3 equations that are all decoupled, which completely misses the point of the interactions between l, M & h in the original system. $\endgroup$ – boschbird Sep 24 '15 at 5:11
  • $\begingroup$ 2) Asymptotics: I can't see the justification for this after 4 years of university mathematics training. $\endgroup$ – boschbird Sep 24 '15 at 5:19
  • $\begingroup$ Thankyou for your response though. $\endgroup$ – boschbird Sep 24 '15 at 5:32
  • 1
    $\begingroup$ @ boschbird: de nada. a) "serious flaw" is a strong expression the use of which you didn't justify. You might know already that simplifying problems can have strong impact on the "physical" content, still it is a common method to come to first results. b) the "justification" is formal, please read again my first sentence c) my linearization for M'[t] == sM (a1*l[t]/B) is correct, you missed l[t] in your comment (this is a "flaw" ;-) $\endgroup$ – Dr. Wolfgang Hintze Sep 24 '15 at 6:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.