Model Description
In the 1830's the Australian Opossum (which is different from the opossums found in the Americas) was introduced to New Zealand for a planned fur trade. However, as is so often the case with this type of species introduction, in New Zealand the opossums had no natural predators, and their numbers increased rapidly. Today they are found throughout New Zealand and are, in addition, a reservoir for bovine tuberculosis—which they may transfer to livestock.[1]
Our model considers the population of opossums, \(P\), which may be broken into a subpopulation of susceptible opossums, \(S\), and those who are infected, \(I\) (thus, \(P = S + I\)). We assume, following [2], that there are no opossums who recover from T.B.
ODE Model
Following [1] and [2], we consider a simple model for the population \(P\) of opossums and the number of those, \(I\), which have contracted tuberculosis, \[\begin{aligned} P' &= (a - b)P - \alpha I\\ I' &= \beta I(P - I) - (\alpha + b) I. \end{aligned}\] In the first of these \(a\) and \(b\) are the natural birth and death rates, respectively, for opossums (we expect \(a > b\), so that in the absence of other effects the population grows exponentially). The parameter \(\alpha\) is the disease-induced death rate, and \(\beta\) is a parameter giving the infection rate; the term \(\beta I(P - I)\) thus represents the infection of animals, which is proportional to contacts between diseased and non-diseased opossums.
We can recast this by rewriting it in terms of the susceptible population by substituting \(P = S + I\) in these equations, to obtain the alternate system \[\begin{aligned} S' &= (a - b)S + aI - \beta I S\\ I' &= -(\alpha + b)I + \beta I S. \end{aligned}\] Non-dimensionalizing by taking \(S = \frac{a-b}{\beta}\, x\), \(I = \frac{a-b}{\beta}\, y\) and \(t = \frac1{a-b}\,s\) we obtain the simpler system \[\begin{aligned} x' &= x + r y - x y \\ y' &= -c y + x y, \end{aligned}\] where \(r = \frac{a}{a-b}\) and \(c = \frac{\alpha+b}{a-b}\). From [1] we expect \(a \approx 1.8 b\), so that \(r \approx \frac94\).
There are two equilibrum solutions for this system, \((x,y) = (0,0)\) and \((x,y) = (c, 4c/(4c -9))\). Linearizing around these by taking \((x,y) = (x_{eq},y_{eq}) + (u,v)\), we have the linear systems \[ \begin{aligned} u' &= u + r v\\ v' &= -c v \end{aligned} \qquad\mbox{and}\qquad \begin{aligned} u' &= -\frac{9}{4c-9}\, u - \frac{4c-9}{9}\, v\\ v' &= \frac{4c}{4c-9}\,u \end{aligned}.\]
In this demonstrations we show the eigenvalues and vectors and graph phase portraits for these. If \(c = 3\) we have the relatively nice eigenvalues and vectors \(\lambda = -\frac32 \pm i\frac{\sqrt{3}}2\), \(\mathbf{v} = \begin{pmatrix} -3\pm i\sqrt3\\ 8\end{pmatrix}\) at the non-zero equilibrium point (and \(\lambda_{1,2} = 1, -3\); \(\mathbf{v}_{1,2} = \begin{pmatrix} 1\\ 0\end{pmatrix}, \begin{pmatrix} -9\\ 16\end{pmatrix}\) at \((0,0)\)).
Matlab Demos
Our demonstrations here show the solutions near each equilibrium solution, and those solutions in the full phase plane.
- Aussie_Opossum.m:
A demonstration that graphs the eigenvectors (if real), direction
field and solution trajectories near either equilibrium point.
Notable configuration parameters include
c_value, the value of \(c\) to consider;show_point, which is set to 0 to show the phase portrait near \((0,0)\) and to 1 to show the non-zero equilibrium point; andinclude_pauses, which determines whether pauses are included between graphs. Note: also requires the files plot_directionfield.m, plot_eigenvectors.m, and plot_solutiontraj.m; and arrow.m (downloads as a zip file with the Matlab file and license). [show figure]
- Aussie_Opossum_Phase_Plane.m:
A demonstration that graphs the phase plane for the system. The
solutions to the linearized system are shown near their respective
equilibrium points, along with trajectories generated numerically
for the full system. Notable configuration parameters include
c_value, the value of \(c\) to consider;show_direction_field, which determines if the direction field is shown as well, andinclude_pause, which indicates whether to pause after graphing the direction field and linear trajectories. Note: also requires the file plot_localtraj.m. [show figure]
Looking at the Model
Some questions that may be worth considering:
- What do the phase portraits near the equilibrium solutions tell us about the behavior of the system for the full phase plane?
- How do we determine the direction of the spiral trajectories that occur when there are complex eigenvalues?
- When we look at the equilibrium solution \((0,0)\), what happens when \(u < 0\) or \(v < 0\)? Is this different at the non-zero equilibrium point?
References
- Wake, G (Winter 1995). "Percy Possum Plunders." CODEE Newsletter. The Community of Ordinary Differential Equations Educators, http://www.codee.org/.
- Louie, K., Roberts, M.G., and G.C. Wake (1993). "Thresholds and Stability Analysis of Models for the Spatial Spread of a Fatal Disease." IMA Journal for Mathematics Applied to Medicine and Biology. 10:207-226.
