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

Question

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?

Answer 1

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] $$