Stability Analysis
Qualitative analysis of dynamical systems involves methods to understand the behavior of the systems without solving the equations. One of the most important tools for mathematicians is stability analysis: we identify steady states or stable orbits, and we investigate whether they are stable. By stable, we mean that the system would tend to return to the steady state after a small perturbation. To examine stability, we look at the dynamical behavior close to an equilibrium. Up close, the dynamics are approximately linear, so we examine the behavior of a linearized system. In the case of developing threshold conditions for malaria persistence, we investigate the stability of the disease free equilibrium.
The origins of \(R_0\) are in human demography, where threshold criteria could be computed as the lead eigenvalue of a life history matrix. The primary function of \(R_0\) was to develop a biologically meaningful formula that could also serve as a threshold. Alfred J Lotka, who developed and named \(R_0,\) also had a long standing interest in malaria [1,2], including a thorough stability analysis of Ross’s equations. Macdonald does not cite Lotka in his papers, but he does give credit to Lotka in his 1957 book [3]. In telling the story, Ross developed the first models, and Macdonald’s synthesis was transformative, but it was Lotka who should get credit for putting the biomathematics of malaria on solid ground.
Here, we take a tour of methods for computing threshold conditions for this basic malaria model:
\[\begin{array}{rl} \frac{dX}{dt} = F_X(Y) &= bfq\frac{Y}{H}(H-X) - r X \\ \frac{dY}{dt} = F_Y(X) &= cfq \frac{X}{H} (M-Y) - g Y \end{array}\]
By now, it should be familiar:
\(X\) is the population density of infected humans among \(H\) total
\(Y\) is the population density of infected mosquitoes among \(M\) total
The bionomic parameters are explained in Notation
Steady States
We note that the system has two steady states. The disease-free steady state is at \(X=Y=0.\)
To compute the endemic equilibrium, let \(x = X/H\) and \(y=Y/M.\) At the steady state,
\[y = \frac{cfqx}{g + cfqx}\]
and substituting this back, we get
\[\frac{dx}{dt} = bcf^2 q^2 \frac{x}{g+cfqx}(1-x) - rx = 0\]
We let
\[R_0 = bc \frac{M}{H} \frac{f^2 q^2}{gr}\]
A recurring term is the expected number of times a mosquito would get infected in its lifetime, assuming humans are perfectly infectious: \[s=cfq/g\]
Since we’re not at \(x=0\), we can divide by \(x,\) rewrite and simplify
\[R_0 (1-x) = 1+s\] Collecting and solving for \(x\), we get
\[ x = \frac{R_0-1}{R_0 + s}\] substituting this back, we get
\[y = \frac{R_0-1}{R_0} \frac{s}{1+s}\]
The steady states are \(\bar X = x H\) and \(\bar Y = y M.\)
Stability Analysis
To evaluate the stability of these equilbria, we start by computing the matrix of partial derivatives:
\[R = \left[ \begin{array}{rl} \frac{\partial F_X}{dX} & \frac{\partial F_X}{dY} \\ \frac{\partial F_Y}{dX} & \frac{\partial F_Y}{dY} \end{array} \right] = \left[ \begin{array}{rl} - bfq\frac{Y}{H} - r& bfq \frac{H-X}{H} \\ cfq\frac{M-Y}{H} & - cfq \frac{X}{H} - g \end{array} \right]\] The question is whether pertubrations would grow. Effectively, we are looking at the stability of the linearized system at these steady states.
If we evaluate \(R\) at the disease-free equilibrium, \(X_0=Y_0=0.\)
\[\left. R \right|_{X=Y=0} = \left[ \begin{array}{rl} - r& bfq \\ cfq\frac{M}{H} & - g \end{array} \right]\] The eigenvalues of a matrix are negative if and only if the trace is negative and the determinant is positive. Since the trace of \(R\) is negative (using mathematical notation, \(\mbox{Tr}\left(R\right) = - r - g\)), the determinant must be positive. The disease free state is stable (meaning the parasite will not increase) only if: \[\mbox{det}\left(R\right) = rg - bc f^2 q^2 \frac{M}{H} > 0 \] It follows that the disease free state is unstable (and the parasite will increase) only if:
\[R_0 = bc \frac{f^2 q^2}{rg} \frac{M}{H} > 1. \]
Eigenvalues & Eigenvectors
We can find the eigenvalues (vectors \(E\) and scalars \(\ell\) that satisfy \(R \cdot x = \ell x\)), by computing
\[\mbox{det}\left[ \begin{array}{rl} - r-\ell& bfq \\ cfq\frac{M}{H} & - g -\ell \end{array} \right] = \ell^2 + \ell \; \mbox{Tr}\left(R\right) + \mbox{det}\left(R\right) = 0 \] or
\[ \ell = \frac{- \mbox{Tr}\left(R\right) \pm \sqrt{\mbox{Tr}\left(R\right)^2 - 4 \; \mbox{det}\left(R\right)} }{2}\]
If we let \(E_1\) be the eigenvector associated with the larger eignvalue (\(\ell_1\)) and \(E_2\) the other one (associated with eigenvalue \(\ell_2\)), then we can write any initial vector \(X_0, Y_0\) as a linear combination:
\[\left[\begin{array}{r} X_0 \\ Y_0 \end{array} \right] = \xi_1 E_1 + \xi_2 E_2 \]
Then
\[e^{Rt} = \xi_1 \ell_1^t E_1 + \xi_2 \ell_2^t E_2\] We will come back to this idea, because it is the basis for some computational methods that shed light on analysis of invasion or reinvasion of malaria in heterogeneous systems.
Next Generation Matrix
Let the infected classes \(\{X,Y\}\) be denoted by a vector \(x_i\). The next generation matrix is computed by splitting the \(F_X(Y)\) and \(F_Y(X)\) into two components: the rate of appearance of new infections \(\mathcal{F(x_i}\) classes in the model and the rate of transfer of individuals in and out of a class \(\mathcal{V(x_i)}\). Here \(V(X_i) = V^-(x_i) - V^+(x_i)\) where \(V^-(x_i)\) is the rate of transfer out of a class and \(V^+(x_i)\) is the rate of transfer of individuals into the class. Thus, \[ \mathcal{F(x_i)} = \left[\begin{array}{r} bfq\frac{Y}{H}(H-X) \\ cfq \frac{X}{H} (M-Y) \end{array} \right], \quad \mathcal{V(x_i)}= \left[\begin{array}{r} r X \\ g Y \end{array}\right]\] Next we determine the matrices \(F =\frac{\partial F(x_i)}{\partial x_i}|_{\{X=0,Y=0\}}\), and \(V =\frac{\partial V(x_i)}{\partial x_i}|_{\{X=0,Y=0\}}\) such that \[ F = \left[\begin{array}{rl} -bfq\frac{Y}{H} & bfq\frac{H-X}{H}\\ cfq \frac{M-Y}{H} & -cfqX\frac{M}{H} \end{array} \right], \quad V= \left[\begin{array}{rl} r & 0 \\ 0 & g \end{array}\right].\] Evaluating at \(X=0, Y=0\) and computing \(FV^{-1}\) yields \[FV^{-1} = \left[\begin{array}{rl} 0 & \frac{bfq}{r}\\ cfq \frac{M}{gH} & 0 \end{array} \right]\]
We let \[R = (F V^{-1})^2 = \left[\begin{array}{rl} bc\frac{f^2 q^2}{rg} \frac{M}{H} & 0\\ 0 & bc\frac{f^2 q^2}{rg} \frac{M}{H} \end{array} \right]\]
In the \(n^{th}\) parasite generation, the number of offspring would be: \[R^n \cdot \left[\begin{array}{r} X_0 \\ Y_0 \end{array} \right] = (R_0)^n \cdot \left[\begin{array}{r} X_0 \\ Y_0 \end{array} \right] \] The numbrer of offspring is growing iff \(R_0 > 1\)
The matrix \(FV^{-1}\) is referred to as the next generation matrix. The approach, which was developed for structured population models by Diekmann et al. and also by van den Driessche and Watmough works for all sorts of diseases [4,5]. In our case, it is the parasite’s generation that matters, and a parasite generation is defined by one full cycle, from human to human, or from mosquito to mosquito.
The spectral radius is commonly interpreted as a basic reproduction number threshold, and in this case, it is called \(R_0^2,\) because each step in transmission –from human to mosquitoes, and from mosquitoes back to humans – is called a generation, and so \[R_0 = \sqrt{bc\frac{f^2 q^2}{rg} \frac{M}{H}}.\] If, however \(R_0 > 1\) is a threshold criterion, then an equally valid cretrion is \(R_0^2 >1.\) The decision to call it \(R_0^2\) is purely semantic. Like Lotka, the formula