# Analytical stability analysis of ODE model does not support numerical integration results

I am analysing the two-dimensional model:

$$\frac{dS}{dt} = \phi S \big(\frac{S}{\tau} - 1\big) \big(1 - \frac{S}{\kappa}\big)$$

$$\frac{dI}{dt} = \beta \rho S - \mu I - \omega I$$

where all parameters are $$>0$$.

I have identified three equilibria:

$$\hat{S} = 0, \hat{I} = 0$$

$$\hat{S} = \tau, \hat{I} = \frac{\beta \rho \tau}{\mu + \omega}$$

$$\hat{S} = \kappa, \hat{I} = \frac{\beta \rho \kappa}{\mu + \omega}$$

I have analysed the stability of these equilibria by computing the Jacobian:

$$\boldsymbol{J} = \begin{bmatrix} \frac{2 \phi S}{\tau} - \phi - \frac{3 \phi S^2}{\kappa \tau} + \frac{\phi S}{\kappa} & 0 \\ \beta \rho & -\mu - \omega \end{bmatrix}$$

and evaluating the Jacobian at the equilibria.

For the second equilibria, I have found that the Jacobian is:

$$\boldsymbol{J} = \begin{bmatrix} \phi - \frac{2 \phi \tau}{\kappa} & 0 \\ \beta \rho & -\mu - \omega \end{bmatrix}$$

Because this matrix is lower-triangular, it has eigenvalues at $$\lambda_{1} = \phi(1 - \frac{2\tau}{\kappa})$$ and $$\lambda_{2} = -\mu - \omega$$.

The second eigenvalue is always negative given $$\mu, \omega > 0$$. The first eigenvalue will be negative if:

$$\phi(1 - \frac{2\tau}{\kappa}) < 0$$

which occurs if $$\tau > \kappa/2$$.

I have tried to verify these results using numerical simulation, but I am finding that the second equilibria is always unstable. My first question is, therefore, is my math correct?

For a numerical example in R:

require(deSolve)

m_twod <- function(t,y,p)
{
with(as.list(c(y,p)), {

dS_dt <- phi*S*(S/tau - 1)*(1 - S/kappa)
dI_dt <- beta*rho*S - mu*I - omega*I

return(list(c(
dS_dt, dI_dt
)
))
})
}

# run the model
run_twod <- as.data.frame(
ode(
func = m_twod,
y = c(S=111, I = 0),
parms = c(phi=0.1, tau=110, kappa = 200,
beta = 1/40, rho=12*110*0.75*(1-0.12), mu = 0.3, omega=0.2),
times = seq(0, 52*5, 0.01),
method = "ode45"
)
)



I ensure that in this model, $$\tau > \kappa/2$$, and start the model near to the $$\hat{S} = \tau$$ equilibrium point. However, the model still goes to the $$\hat{S} = \kappa, \hat{I} = \beta \rho \kappa/(\mu+\omega)$$ equilibria.

Could anyone shed some light on my analysis problems?

According to my calculations from

$$\left\{ \begin{array}{rcl} \phi S(t) \left(1-\frac{S(t)}{\kappa }\right) \left(\frac{S(t)}{\tau }-1\right) & = & 0\\ \beta \rho S(t)-(\mu +\omega ) I(t) & = & 0\\ \end{array} \right.$$

we have the three equilibrium points

$$\left[ \begin{array}{ccc} n & S & I\\ 1 & 0 & 0 \\ 2 & \kappa & \frac{\beta \kappa \rho }{\mu +\omega } \\ 3 & \tau & \frac{\beta \rho \tau }{\mu +\omega } \\ \end{array} \right]$$

The Jacobian gives

$$J = \left( \begin{array}{cc} -\frac{\phi \left(3 S^2-2 (\kappa +\tau ) S+\kappa \tau \right)}{\kappa \tau } & 0 \\ \beta \rho & -(\mu +\omega) \\ \end{array} \right)$$

and at each equilibrium point the eigenvalues

$$\left[ \begin{array}{ccc} n & & \\ 1& -(\mu +\omega) & -\phi \\ 2& \phi -\frac{\kappa \phi }{\tau } & -(\mu +\omega) \\ 3& -(\mu +\omega) & \phi -\frac{\phi \tau }{\kappa } \\ \end{array} \right]$$

• Thanks! I realised that the first partial derivative in the Jacobian was wrong. Commented Mar 24, 2020 at 15:07