Chapter 4 Equilibrium Stability Analysis
4.1 Introduction to the Chapter
In Chapter 3.4, we briefly introduced equilibrium and \(R_0\) (Definition 1.3) for examining stability in equilibrium (Definition 3.2). Chapter 4 will expand on the topic of equilibrium through linear stability analysis.
Definition 4.1 (Linear Stability Analysis) A method providing the local stability of each equilibrium, as we approximate the solutions near the equilibrium. The approach to linear stability analysis involves determining the eigenvalues of the Jacobian matrix at the given equilibrium. We can linearize the system of equations around the equilibria with the Jacobian matrix and find the characteristic polynomial whose roots are the eigenvalues. We will show that all the eigenvalues have negative real parts to maintain the local stability of the equilibrium. See Chapter 1.2 for the recommended linear algebra and calculus background.
This chapter will explore both the direct calculation of eigenvalues, as well as some alternative / indirect methods of showing stability.
4.2 Linear Stability Analysis for the SIR Model
The stability of an equilibrium in a compartmental model (Definition 2.3) is determined by the signage of the eigenvalues (\(\lambda\)) of the Jacobian Matrix (J). In this section, we will solve for the eigenvalues of the SIR model at the Disease Free Equilibrium (DFE) and the Endemic Equilibrium (EE).
4.2.1 Direct Calculation
The SIR model has 3 compartments, and thus, 3 differential equations. The equation below represents the balance of each compartment, totaling to 1. We can rearrange this to substitute R out, as shown below.
\[S+I+R=1 \Rightarrow R=1-S-I\] Now, we have our 3 system of equations, (\(\dot{S}, \dot{I}, \dot{R}\)), in which \(\dot{S}, \dot{I}\) represent Functions F and G respectively. By using the equation above we are able substitute R, and reduce our system from 3 to two equations, (\(\dot{S}, \dot{I}\)) \[ \begin{align*} \dot{S} &= F(S,I) = -\beta SI + \alpha R \\ \dot{I} &= G(S,I) = \beta SI - \gamma I \\ \dot{R} &= \gamma I - \alpha R \\ \end{align*}\]
Our substitution results in these 2 equations for \(\dot{S}, \dot{I}\):
\[\begin{align} \dot{S} &= -\beta SI + \alpha(1-S-I) \\ \dot{I} &= \beta SI - \gamma I \end{align}\]
The Jacobian Matrix (J) represent a matrix of the partial derivatives of the system of Functions F and G. As stated before, this is represented by \(\dot{S}, \dot{I}\)
\[J = \begin{bmatrix} \frac{dF}{dS} & \frac{dF}{dI} \\ \frac{dG}{dS} & \frac{dG}{dI}\end{bmatrix}\] By substituting our system of equations, where \(S^*\) and \(I^*\) are the values of S and I at our specific equilibrium point, the Jacobian Matrix is as follows:
\[J = \begin{bmatrix}-\beta I^* - \alpha & -\beta S^*-\alpha \\ \beta I^* & \beta S^* - \gamma\end{bmatrix}\]
4.2.1.1 At Disease-Free Equilibrium
The Disease-Free Equilibrium (DFE) (Definition 3.3) is the state in which the disease has been eradicated. Therefore, at DFE, \(I^* = 0\) and \(S^* = 1\). It is important to note that this equation must retain the balance from our earlier equation (\(S+I+R = 0\)), so \(R^* = O\).
We can substitute our values for \(I^*\) and \(S^*\) at the DFE into our Jacobian Matrix before, so we get:
\[J_{DFE} = \begin{bmatrix}-\alpha & -\beta-\alpha \\ 0 & \beta-\gamma\end{bmatrix}\]
Next, we solve for the eigenvalues of the Jacobian Matrix by the formula \(det(J-\lambda I_n) = 0\), where \(det\) is the determinant, \(J\) is the Jacobian Matrix, \(\lambda\) is the eigenvalue, and \(I_n\) is the identity matrix with \(n\) dimensions.
Adding in the \(J\) we solved for above and the identity matrix \(I_2 = \begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix}\), we get:
\[ 0 = det(\begin{bmatrix}-\alpha & -\beta-\alpha \\ 0 & \beta-\gamma\end{bmatrix}-\lambda\begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix}) \] We will now simplify the contents of the determinant using matrix operations.
\[ \begin{align*} \ 0 &= det(\begin{bmatrix}-\alpha & -\beta-\alpha \\ 0 & \beta-\gamma\end{bmatrix}+\begin{bmatrix}-\lambda & 0 \\ 0 & -\lambda\end{bmatrix}) \\ \ 0 &= det(\begin{bmatrix}-\alpha-\lambda & -\beta-\alpha \\ 0 & \beta-\gamma-\lambda\end{bmatrix}) \\ \end{align*} \] Evaluate the determinant of the matrix by the rule \(det \left(\begin{bmatrix} A & B \\ C & D\end{bmatrix} \right) = AD-BC\). Applying this rule, we get:
\[ 0 = -(\alpha+\lambda)(\beta-\gamma-\lambda) \] Then, solve for the two eigenvalues.
\[ \begin{align*} \ 0 =-(\alpha+\lambda_1) &\longrightarrow \lambda_1 = -\alpha\\ \ 0 =(\beta-\gamma-\lambda_2) &\longrightarrow \lambda_2 = \beta-\gamma \\ \end{align*} \]
Therefore, the eigenvalues are \(\lambda_1 = -\alpha\) and \(\lambda_2 = \beta-\gamma\).
Definition 4.2 (Equilibrium Stability by Eigenvalues)
- The equilibrium is stable if all the eigenvalues have negative real parts.
- The equilibrium is unstable is any eigenvalue is positive.
- The equilibrium stability is inconclusive if some eigenvalues are zero and the rest are negative.
Stability Analysis of eigenvalues at DFE using Definition 4.2: For DFE to be stable, we will show both the eigenvalues \(\lambda_1\) and \(\lambda_2\) have negative real parts.
- Showing \(\lambda_1\) has negative real parts: \[(\lambda_1 = -\alpha) < 0\]
- Showing \(\lambda_2\) has negative real parts: \[ \begin{align*} \ (\lambda_2 = \beta - \gamma) &< 0 \\ \ \beta &< \gamma \\ \ \frac{\beta}{\gamma} &< 1 \\ \ R_0 &< 1 \\ \end{align*} \]
Therefore, we have shown that the eigenvalues have negative real parts and also obtain that the Disease Free Equilibrium is stable when \(R_0 = \frac{\beta}{\gamma} < 1\)
4.2.1.2 At Endemic Equilibrium
The Endemic Equilibrium (Definition 3.4) represents the state in which the disease persists. Recall from Chapter 3.4 that at EE, we solve the system of equations (\(\dot{S}, \dot{I}\)) where \(R = 1-S-I\) and obtain (\(S^*\), \(I^*\)).
\[ \begin{align*} \ \dot{S} &= -\beta SI + \alpha(1-S-I) &\Longrightarrow S^* &=\frac{\gamma}{\beta} \\ \ \dot{I} &= \beta SI - \gamma I &\Longrightarrow I^* &= \frac{\alpha (1 - \frac{\gamma}{\beta})}{\gamma + \alpha} \\ \end{align*} \]
Note that since \(S^*=\frac{\gamma}{\beta}=\frac{1}{R_0}\), we have \(R_0>1\). Thus, if EE exists, the DFE is unstable.
We substitute our value for \(S^*\) into the Jacobian Matrix we defined at the start of Chapter 4.2.1. For simplicity, we will keep \(I^*\) as is for now.
\[ \begin{align*} \ J &= \begin{bmatrix}-\beta I^* - \alpha & -\beta S^*-\alpha \\ \beta I^* & \beta S^* - \gamma\end{bmatrix} \\ \ J_{EE} &= \begin{bmatrix}-\beta I^* - \alpha & -\gamma-\alpha \\ \beta I^* & 0\end{bmatrix} \\ \end{align*} \]
To solve for the eigenvalues of the Jacobian Matrix, we follow the same approach used earlier for the DFE. Use the formula \(det(J-\lambda I_n) = det(\lambda I_n -J) = 0\), and further simplify by matrix operations and determinant rules. We obtain the following:
\[ \begin{align*} \ 0 &= det \left(\begin{bmatrix}-\beta I^* - \alpha -\lambda & -\gamma-\alpha \\ \beta I & -\lambda\end{bmatrix} \right) \\ \ 0 &= -(\beta I^* + \alpha + \lambda) (-\lambda) + (\beta I^*) (\gamma+\alpha) \\ \ 0 &= \lambda^2 + (\beta I^* + \alpha)\lambda + \beta I^*(\gamma+\alpha) \\ \end{align*} \] Finally we substitute \(I^*\) into the equation and get: \[ \begin{align*} \ 0 &= \lambda^2 + \left[ \beta \left(\frac{\alpha (1 - \frac{\gamma}{\beta})}{\gamma + \alpha} \right) + \alpha \right] \lambda + \beta \left( \frac{\alpha (1 - \frac{\gamma}{\beta})}{\gamma + \alpha} \right) (\gamma+\alpha) \\ \ 0 &= \lambda^2 + \alpha \left[ \frac{1 - \frac{1}{R_0}}{R_0 + \frac{\alpha}{\beta}} + 1 \right]\lambda + \alpha\beta \left(1-\frac{1}{R_0} \right) \\ \end{align*} \]
The equation above is the characteristic equation, a polynomial equation produced from the Jacobian Matrix. This characteristic equation can be solved to find the eigenvalues.
In this case, because the characteristic equation is a quadratic equation, we can apply the quadratic formula to solve for the eigenvalues.
Recall the quadratic formula: \[x = \frac{-b \pm \sqrt{b^2 - 4ac}} {2a}\]
We have \(a = 1\), \(b = \alpha \left[ \frac{1 - \frac{1}{R_0}}{R_0 + \frac{\alpha}{\beta}} + 1 \right]\), and \(c=\alpha\beta \left(1-\frac{1}{R_0} \right)\).
Substituting a, b, and c into quadratic formula:
\[\lambda = \frac{-\alpha \left[ \frac{1 - \frac{1}{R_0}}{R_0 + \frac{\alpha}{\beta}} + 1 \right] \pm \sqrt{(\beta I^* + \alpha)^2 - 4\beta I(\gamma+\alpha)}}{2}\]
The general solution to the eigenvalues is represented by \(\lambda = \frac{-b}{2}\) \(\pm\) Complex Number, where \(\frac{-b}{2}\) is the real part. Therefore we get the two eigenvalues of the form \(\lambda_1 = v + wi\) and \(\lambda_2 = v-wi\), where \(Re(\lambda)=v\) represents the real part, \(Im(\lambda)=w\) represents the complex part, and \(i = \sqrt{-1}\) is the imaginary number. The complex part only exists if the sign of the discriminant is negative.
Stability Analysis of eigenvalues at EE using Definition 4.2: For EFE to be stable, we will show both the eigenvalues \(\lambda_1\) and \(\lambda_2\) have negative real parts. Since eigenvalues of EE may potentially have a complex part, we must consider both the real and complex parts when interpreting stability.
We know that the real part is negative from the signage of \(b = -\alpha \left[ \frac{1 - \frac{1}{R_0}}{R_0 + \frac{\alpha}{\beta}} + 1 \right]\). By Definition 4.2 we know this indicates that the EE is stable.
To further our understanding, we can consider the case where the complex part exists. If \(c > 0\), then \(0 < b^2 - 4c < b^2\). Recall that \(a = 1\). Then:
\[ \begin{align*} \lambda_1 &= \frac{{-b} - \sqrt{b^2 - 4c}}{2} < 0 \\ \lambda_2 &=\frac{{-b} + \sqrt{b^2 - 4c}}{2} < \left( \frac{{-b} + \sqrt{b^2}}{2}=0 \right) \\ \end{align*} \]
- Therefore we have shown that the eigenvalues are negative when the complex part is present or absent, so the Endemic Equilibrium is stable when c>0.
4.2.2 A Note on Phase Diagrams
A phase diagram is a visual representation of the conditions that can occur at equilibrium.
The phase diagram for the SIR model consists of a simplex where the diagonal line represents the boundary of no recovered individuals, and the origin represents the point at which all individuals are recovered. The space within the simplex is when there are some recovered individuals. The equilibrium exists somewhere within this simplex.
Different phase diagrams can exist depending on the conditions:
- When \(R_0 < 1\), a single equilibrium exists: DFE. The phase diagram will have a stable node, in which trajectory of the solution goes towards the equilibrium.
- When \(R_0 > 1\), two equilibrium exist: DFE and EE. DFE is always an unstable node and EE is either a stable node or a stable spiral The trajectory of the solution goes from the unstable DFE to the stable EE. Whether EE is a spiral depends on whether a complex solution exists for the eigenvalues. If the discriminant is negative, a complex value is present and a stable spiral will occur. If the discriminant is nonnegative, there is no complex number and a stable node will occur.
4.2.3 Routh-Hurwitz Criteria
For stability, it must be shown that the real part of all the eigenvalues is negative (Definition 4.2).
Definition 4.3 (Routh-Hurwitz Criteria) A method to check if the eigenvalues are negative and the system is stable without needing to solve for the roots. The criteria can determine if the eigenvalues are negative by looking at the coefficients. Depending on the degree of the polynomial, the criteria varies.
We will consider stability analysis for different characteristic polynomials using the Routh-Hurwitz Criteria (Definition 4.3).
4.2.3.1 Quadratic Characteristic Equation
For a quadratic characteristic equation, \(c_2\lambda^2 + c_1\lambda + c_0=0\), the coefficients must all be positive. As long as \(c_2\), \(c_1\), and \(c_0\) are positive the eigenvalues will be negative.
The quadratic characteristic equation can be used when there are 3 compartments. The SIR model, with the susceptible, infected, and recovered compartment, can be used with the quadratic characteristic equation.
4.2.3.2 Cubic Characteristic Equation
For a cubic characteristic equation, \(c_3\lambda^3+c_2\lambda^2 + c_1\lambda + c_0=0\), the criteria are slightly more complicated. Here the conditions that must be met are :
\[ c_i>0 \\ c_1c_2>c_0c_3 \\ \] And since \(c_3 = 1\), the conditions are therefore: \[ c_i>0 \\ c_1c_2>c_0 \\ \] As long as these conditions are met, it can be found that a cubic characteristic equation has only negative real parts and the system is therefore stable.
A cubic characteristic equation can be used for models with 4 compartments such as the SEIR (Susceptible, Exposed, Infected, Recovered) model, SAIR (Susceptible, Asymptomatic, Infected, Recovered) model, and the SIVR(Susceptible, Infected, Vaccinated, Recovered) model.
4.2.3.3 Quartic Characteristic Equation
For a quartic characteristic equation, \(c_4\lambda^4+ c_3\lambda^3+c_2\lambda^2 + c_1\lambda + c_0=0\), the criteria are further complicated. Here the conditions that must be met are:
\[ c_i>0 \\ c_2c_3>c_1c_4 \\ c_1c_2c_3 - c_4c_1^2 > c_0c_3^2 \] As long as these conditions are met, it can be found that a quartic characteristic equation has only negative real parts and the system is therefore stable.
4.2.4 Descartes’ Rule of Signs
Definition 4.4 (Descartes' Rule of Signs) A method to determine the number of potential real positive and/or negative roots.
We will apply Descartes’ Rule of Signs (Definition 4.4) for stability analysis to the following characteristuc equation.
\[c_2\lambda^2 + c_1\lambda + c_0 = 0\]
Using this general characterisitic equation, as seen when finding the eigenvalues of the Jacobian Matrix, we can understand how stable the equilibrium is in the system of differential equations. For stability to occur, all real roots and the real part of complex roots must be negative (Definition 4.2).
4.2.4.1 Case 1: Positive Roots with No Sign Changes
\[(+)(+)(+)\]
0 sign changes from term to term
All terms are positive, denoted by each coefficient (\(c_0, c_1, c_2,\)) being positive.
- Result: There are either 0 positive roots
4.2.4.2 Case 2: Negative Roots with Multiple Sign Changes
\[\begin{align*} c_2(-\lambda)^2+c_1(-\lambda)+c_0=0 \\ c_2\lambda^2-c_1\lambda+c_0=0 \end{align*}\]
- 2 sign changes from term to term
- First term (\(c_2\lambda^2\)) and last term (\(c_0\)) are positive, while middle term (\(c_1\lambda\)) is negative
- From first to second term, \(+\) to \(-\), first sign change
- From second to third term, \(-\) to \(+\), second sign change
- Number of real negative roots = Number of sign changes
- Result: There are either 2 negative roots or 0 negative roots
4.2.4.3 Cases with One Sign Change:
Example:
\((-) \rightarrow (+) (+)\)\(\cdots\)\(\lambda \rightarrow\) 1 Positive Root
- 1 Sign Change and Positive Coefficients = 1 Real Positive Root
\((-) (-) \rightarrow (+)\)\(\cdots\)\(-\lambda \rightarrow\) 1 Negative Root
- 1 Sign Change and Negative Coefficients = 1 Real Negative Root
4.3 Linear Stability Analysis for the SI Model
An additional example is provided for the stability analysis of the SI model. \(R_0\) is used to analyze the stability of this simpler model.
\[ \begin{align*} \ \dot{S} &= -\beta SI + \alpha I =(-\beta S \gamma)(1-S)\\ \ \dot{I} &= \beta SI - \gamma I\\ \end{align*} \]
Solving the Jacobian:
\[
det(J) = -\beta(1-S) - (-\beta S + \gamma) = 2\beta S - \beta - \gamma
\]
At DFE, \(I^*=0\), \(S^*=1\), so:
\[
det(J_{DFE}) = 2\beta (1) - \beta - \gamma = \beta - \gamma
\]
\[ \begin{align*} \ \beta - \gamma &< 0 \\ \ \frac{\beta}{\gamma} &< 1 \\ \ R_0 &< 1 \\ \end{align*} \]
At EE, \(S=\frac{\gamma}{\beta}\), so: \[ det(J_{EE}) = 2\beta \left( \frac{\gamma}{\beta} \right) - \beta - \gamma = \gamma - \beta \]
\[ \begin{align*} \ \gamma -\beta &< 0 \\ \ \frac{\beta}{\gamma} &> 1 \\ \ R_0 &> 1 \\ \end{align*} \]