6
$\begingroup$

I have the following coupled non-linear ODE system, which describes a biological system:

$$ \begin{cases} \dfrac{dp}{dt} = -\gamma p f,\\ \\ \dfrac{df}{dt} = \gamma p f,\\ \\ \dfrac{dT}{dt} = \left( 1 - \dfrac{k}{T} \right) \left\{ b (1 - T) \dfrac{r^n}{f^n + r^n} - m \dfrac{f^n}{f^n + r^n} \right\}, \end{cases} $$ where the parameters are: $\gamma = 0.5432$, $k = 0.0026$, $b = 0.0885$, $n = 14.9832$, $r = 0.8265$, and $m = 0.9780$.

I want to investigate the stability of the system. To this end, first, one observes that the first two equations imply that $p + f$ must be a constant, let's call it $w$, and in our case, it must be non-negative (for the chosen parameters' values of the system in the above, it's $0.801$). Now, by excluding $p$ from the above system, we can write:

$$ \begin{cases} \dfrac{df}{dt} = \gamma (w - f) f,\\ \\ \dfrac{dT}{dt} = \left( 1 - \dfrac{k}{T} \right) \left\{ b (1 - T) \dfrac{r^n}{f^n + r^n} - m \dfrac{f^n}{f^n + r^n} \right\}. \end{cases} $$

The second system has four equilibria, $(f, T)$:

$$(0, k), (0, 1), (w, k), \left( w, 1- \frac{m}{b} (w / r)^n \right).$$

Now, by calculating the eigenvalues of the Jacobian, it turns out that the first two equilibria are unstable, and the last two are stable equilibria.

My question is as follows: If we simulate our original, i.e., the first system, using Mathematica for some initial conditions:

s = NDSolve[{p'[t] == -gamma p[t] f[t], f'[t] == gamma p[t] f[t], 
T'[t] == b (1 - T[t]) (r^n/(f[t]^n + r^n)) (1 - 0.0026/T[t]) - m (f[t]^n/(f[t]^n + r^n)) (1 - 0.0026/T[t]), 
p[0] == p0, f[0] == 0.001, T[0] == T0}/.{gamma -> 0.5432, b -> 0.0885, m -> 0.9780, r -> 0.8265, n -> 14.9832, p0 -> 0.8, T0 -> 0.8}, {p, f, T}, {t, 0, 720}];

Plot[Evaluate[{f[t], T[t]}/.s], {t, 0, 720}, PlotStyle -> {Blue, Red}, AxesOrigin -> {-1, 0}, PlotRange -> All]

we obtain:

enter image description here

From the above plot it's clear that solutions are converging to $(w, k)$ in the long run, which here $w$ is around $0.8$. However, it seems that the system never converges to $\left( w, 1- \frac{m}{b} (w / r)^n \right)$, for any initial conditions. What is the reason behind that? Can a system never converge to one of its stable equilibria?

EDIT

The reason for transforming the original system: If one wants to do the stability analysis for the original system, one runs into a problem. It has two equilibria: $(f \to 0, T \to k)$ and $(f \to 0, T \to 1)$, regardless of the value of $p$. Since $p$ can be anything, it isn't clear how one can calculate the Jacobian. However, by transforming the system, we can eliminate this issue.

$\endgroup$
7
  • $\begingroup$ $f'=\gamma w-\gamma f^2$ can be solved explicitly with separation of variables and its friends (you need special attention for $\gamma w=0\lor f^2(t_0)= w$). $\endgroup$ Commented Mar 20 at 9:48
  • 1
    $\begingroup$ @Anovice What happens if your initial conditions are very close to $(w,1-(m/b)(w/r)^n)$? $\endgroup$
    – user
    Commented Mar 20 at 10:26
  • 1
    $\begingroup$ "Can an ODE system never converge to its stable equilibrium in the long run?" In short, YES. E.g. the harmonic oscillator oscillates forever and never converges to equilibrium. $\endgroup$
    – Quillo
    Commented Mar 20 at 12:54
  • 1
    $\begingroup$ There is no issue with setting $w=p+f=$ constant. Then $w$ appears simply as a parameter in your new system. At the end you must remember to set $w=p(0)+f(0)$, and that's all. Your last equilibrium point is stable and solutions asymptotic to it can be found simply by setting $T(0)=-5$ in your code. For some very specific choice of parameters, the Jacobian may become singular indicating that the linear term in the approximation near the equilibria vanishes... but that isn't the case with your given numbers $\endgroup$
    – Sal
    Commented Mar 20 at 19:25
  • 1
    $\begingroup$ @Anovice no problem :) $\endgroup$
    – Sal
    Commented Mar 20 at 23:14

1 Answer 1

4
$\begingroup$

The direct answer to your question is "no". At least if by "stable" you mean all eigenvalues of the Jacobian at $x^*$ have negative real part. (See, e.g., Section 2.9 Theorems 1 and 2 from Differential Equations and Dynamical Systems, L. Perko 3rd edition).

How then do we square this with this particular example? The issue comes when you simplify the system. You are correct in saying $w=f+p$ is constant, but you have to understand what this means: This means that as time evolves, $w$ is constant for a given trajectory it does not mean that all trajectories have the same constant $w$. Thus to successfully transform the system we would need to include

$$\frac{dw}{dt} = 0$$

Note that from the original system it is clear that the last two fixed points you give cannot be fixed points of the original system since the first two equations imply right away that $f=0$ or $p=0$ for any fixed point.

Edit: I decided to add some information about the system since this is an interesting case. We have that either $p=0$ or $f=0$. Start with $f=0$. Then we have a fixed point for any value of $p$ and $T=k$ or $T=1$. If instead $p=0$ then either $T=k$ (and $f$ is any value) or $T=1-\frac{mf^n}{br^n}$. Thus we get four curves which fully describe all the fixed points of the system. We can still analyze these fixed points as we would anyway. The issue is that these fixed points are what is called "nonhyperbolic"--you will see that the Jacobian of each has a zero eigenvalue. As the book I referenced puts it, "The question as to whether a nonhyperbolic equilibrium point is stable, asymptotically stable, or unstable is a delicate question." The reason is that linearization "doesn't work" in the nonhyperbolic case.

$\endgroup$
6
  • 1
    $\begingroup$ Thank you for your additional information. The fact remains that the transformed system you give is unfortunately not equivalent to the original system, so analyzing the new system will not give any, or at least complete, information about the original system. You will need to look into the theory of so-called non-isolated fixed points and study the original system. $\endgroup$ Commented Mar 20 at 15:10
  • 1
    $\begingroup$ @Anovice, I added some information to my answer which should help to clarify why the issues you are touching on are unavoidable. $\endgroup$ Commented Mar 20 at 15:27
  • 1
    $\begingroup$ I think this is not quite right, or at least not relevant. You say the last two fixed points in OP cannot be so- but I think you're mistaken: if $p=0$ then $f=w$. Also, those fixed points have nonvanishing linear coefficient when expanding the EOMs about them. Check it and see- just a little differentiating ;) $\endgroup$
    – Sal
    Commented Mar 20 at 21:39
  • 1
    $\begingroup$ Thanks, you are correct that those two points are indeed fixed points of the original system within a certain way of viewing it. Still it's not clear since $w$ is being treated as constant. Still OP has incorrectly categorized the fixed points. There are no isolated fixed points of the system. For instance, at $(p,f,T) = (0,f,k)$, I get that $\frac{\partial}{\partial f}\frac{dp}{dt}=\frac{\partial}{\partial f}\frac{df}{dt} = \frac{\partial}{\partial f}\frac{dT}{dt} = 0$. If you mean to discuss the "simplified system" again that is a waste of time since it is not actually equivalent. $\endgroup$ Commented Mar 21 at 17:37
  • 1
    $\begingroup$ The 'certain way of viewing it' is that all the time derivatives vanish, is there any other? The transformed system is equivalent to the original if $w$ is chosen consistently with the initial conditions. If you study the original system, and find a zero eigenvalue, you say the point is indeterminate (the 'delicate question' you say). To classify it requires looking at higher order terms or using a conserved quantity, which OP has already done (also they have constructed the bijection between solutions of the original and transformed systems). Quite a lot for what you call a 'waste of time' $\endgroup$
    – Sal
    Commented Mar 21 at 22:50

You must log in to answer this question.