# To solve or not to solve, that is the problem.

The senior school Mathematics syllabus is often restricted to the study of single variable differential equations of the first order. Unfortunately most real life examples do not follow such types of relations. In addition, very few differential equations in real life have exact solutions that can be expressed in finite terms (Jordan & Smith, 2007, p. 2). Even if the solution can be found exactly it may be far too difficult to be clearly articulated such as those that form an infinite series. In either case, these real life problems are well beyond the scope of the secondary student to solve.

Does this mean that many of the exciting relationships and models found in the real world cannot be studied by the secondary student? What if we were not interested in the solution? Wait a minute ... What did you say? A mathematician not interested in the solution! Surely that is not right! What if the behaviour of the solution was just as important as the solution itself? What device can be so powerful? Enter the phase plane--a geometrical device.

To understand how the phase plane works, we will first consider the predator-prey model defined by Alfred Lotka in 1920 and Vito Volterra in 1926 called the Lotka-Volterra System (as cited in Murray, 2002, p. 79).

Predator-prey model

The predator-prey model is described using two coupled differential equations

dx/dt = [alpha]x - [beta]yx

dy/dt = [delta]xy - [gamma]y (1)

where:

* x is the number of prey;

* y is the number of predators;

* dt/dx and dy/dt are their respective population growth rates; and

* [alpha], [beta], [gamma], [delta] are parameters representing the interaction between the species:

--[alpha]x represents the growth rate of the prey in the absence of predators;

--[beta]xy represents the rate of predation upon the prey;

--[gamma]y represents the loss rate of predators due to death or migration; and

--[delta]xy represents the growth rate of predators (Liu, Zhang, & Chen, 2005).

It is obvious that the solution to these differential equations is beyond the scope of the senior school Mathematics course. As can be seen, the rates of growth for both predators and prey depend on the population of both the predator and prey at any given time. For this reason, these equations are said to be coupled. Neither equation can be examined independently from the other.

To obtain a relation between y and x we remove the time parameter by

dy/dx = (dy/dt)/(dx/dt) = [delta]xy - [gamma]y/[alpha]x - [beta]yx

If we choose an initial population of prey and predator arbitrarily ([x.sub.0], [y.sub.0]), the direction in which the predator and prey populations will move can be determined by the gradient. By using a numerical method it is possible to construct a plot of the phase path or trajectory that this population will follow as it slides in the direction of the gradient. The entire pattern of phase paths is called the phase plane or phase portrait (Jordan & Smith, 2007, p. 6). Various free to use software tools are available on the World Wide Web to generate such phase planes. The Texas Instruments TI-Nspire graphics calculator can also be used to generate phase portraits (TI-Nspire Student Software Guide, p. 267).

According to Mahaffy (2001) the following approximations for the interaction of lynx (the predators) and the snowshoe hare (the prey) in the Canadian forests were determined: [alpha] = 0.4, [beta] = 0.018, [gamma] = 0.8 and [delta] = 0.023. The chosen time scale affects the parameter as they affect associated growth rates. In this example the time scales are in years. The interaction between lynxs and hares were chosen as a result of extensive population records from the Hudson Bay Company from the early 1800s to 1900s. An example of a phase plane for the interaction of lynx and snowshoe hare populations is shown in Figure 1.

Figure 1 was generated using several different starting populations. It can be seen that the nature of interaction between species is cyclic. Also, for a specific starting population at (34, 23), the population of predator and prey is stable.

Students generating such a phase plane may be required to give reasons to explain the cyclic nature of such a population growth. As the prey is consumed by the predators the predator population increases and the prey population decreases. Since the decrease in prey population cannot sustain the increased population of predators the predator population starts to decrease. This then allows the population of prey to increase due to the decreased population of predators, and so the cycle continues.

[FIGURE 1 OMITTED]

Students also possess the necessary skills to determine points where neither species populations are changing, that is points where their respective growth rates are zero. These points are called nullclines where

dx/dt = 0 and dy/dt = 0

From (1) nullclines for the prey are found from setting [alpha]x - [beta]yx = 0 and their geometric representations are defined by y = [alpha]/[beta] and x = 0. For the predator [delta]xy - [gamma]y = 0 and provide nullclines where x = [gamma]/[delta] and y = 0. A steady state occurs at the intersection of these nullclines so that neither population are changing with time. These occur at (0, 0) the trivial case and

([gamma]/[delta], [alpha]/[beta])

the point where both populations are in perfect harmony.

This is a perfect example to introduce students to such pairs of coupled differential equations. By looking at the behaviour of the solution rather than the exact solution secondary students are able to analyse and predict behaviours of systems which are far more complex than those generally studied in the secondary course.

Infectious disease model

We will now look at a more complex system of differential equations. In 1927, McKendrick and Kermack (as cited in Murray, 2005) began work on the first paper that provided the differential equations for a deterministic general epidemic. According to Murray (2005, p. 245), the interaction between susceptible, infected and recovered individuals, called the SIR model, is described by the following differential system of equations.

dS/dt = -[beta]SI + [gamma]R

dI/dt = [beta]SI - [upsilon]I

dR/dt = [upsilon]I - [gamma]R (2)

where:

* S is the number of susceptible persons to the disease;

* I is the number of infected persons;

* R is the number of recovered individuals;

* the parameter P is the rate of infection;

* [gamma] is the rate of re-infection;

* [upsilon] is the rate at which infection is removed.

The epidemic model proposed by McKendrick and Kermack is the above case where[gamma] = 0, where recovered individuals cannot be re-infected (Edelstein-Keshet, 2005).

In this model the total population N is considered constant such that N = S + I + R. Since N is constant, we can eliminate R from equations by letting R = N - S - I, so that the new set of differential equations are

dS/dt = -[beta]SI + [gamma]R (N - S - I)

dI/dt = [beta]SI - [upsilon]I (3)

Mkhatshwa and Mummert (2010) investigated the use of the SIPR model for infectious diseases in the 2003 outbreak of Severe Acute Respiratory Syndrome (SARS). The SIPR model is identical to SIR model however individuals labelled super-spreaders, P, are included with an additional differential equation for

dP/dt

This SIR model does not include such individuals. During the study approximations were found for the SIR parameters, [beta] = 0.01, [Upsilon] = 0.25, and [gamma] = 0.1. Their research showed that the inclusion of super spreaders produced improved modelling results for the SARS epidemic than the classical SIR model. It also demonstrated that without any control measures the outbreak of SARS would have been extreme.

Figure 2 displays a phase portrait of the SIR model as applied to the SARS outbreak with a total of 800 persons. This paper is not arguing that the SIR model best describes the SARS virus, but rather how does the SARS outbreak behave if it follows the SIR model. It is readily seen that no matter what the starting populations for infected and susceptible persons is, the number of infected patients will eventually reach a steady state of 222. This steady state is said to be stable.

It can be shown using the same method as used in the predator-prey model that there are two steady states. A trivial steady state with S = 0 and I = 0 and a second more meaningful steady state at

S = v/[beta] and I = [gamma] (N - [upsilon]/[beta])/[upsilon] + [gamma]

[FIGURE 2 OMITTED]

This second steady state introduces an important condition on the total population size. If N < [upsilon]/[beta], the steady state for the infected population becomes negative (impossible) and hence an infection does not take hold.

Let us now pose a problem to the senior secondary mathematics student. Assume a research grant has been given to a group of medical researchers to investigate methods to reduce the number of infected SARS patients and that it transmits according to the deterministic equations for epidemic growth with the same parameters listed above. Due to funding and staffing levels, the research team are only able to focus their attention on the recovery rate, infection rate, or preventing re-infection. The student is to give advice with justification on which parameter they should focus.

The objective in this task is to reduce I in the long term. Mathematically we are trying to reduce the steady state where

I = [gamma] (N - [upsilon]/[beta])/([upsilon] + [gamma])

However, the laboratory can only focus on changes to one of the parameters [beta], [gamma] and [upsilon]. Reductions in I require that the infection rate [beta] and re-infection [gamma] be both reduced whilst the rate at which infection is removed, [upsilon], must be increased.

In order to compare which parameter causes the greatest reduction in

I = [gamma](N - [upsilon]/[beta]/([upsilon] + [gamma])

we change one parameter whilst keeping the remaining parameters constant. We can then compare the steady state for the number of infected patients to determine which parameter caused the lowest infected population. Since the change is not likely to be linear it is necessary to test for a variety of percentage changes.

Figure 3 presents a graph displaying percentage changes for [beta], [gamma], [upsilon] versus final stable infected population, I. It is evident from this graph that a significant recommendation may be given to the research laboratory. Clearly the research laboratory should focus on reducing the rate of re-infection. This may be counter-intuitive as one would expect the focus to be on minimising infection rates. Hence, this mathematical analysis would prevent precious time wasting research.

[FIGURE 3 OMITTED]

This qualitative analysis of differential equations was founded by Henri Poincare and Ivar Bendixson (cited in Jordan & Smith, 2007, p. 405). These are but a few examples of dynamical systems and the reader is encouraged to explore other such systems. There are many other interesting aspects of dynamical systems just waiting to be investigated such as bifurcations where the changing of certain parameters cause catastrophic changes in behaviour.

Although methods used to find exact solutions to these coupled differential equations may be beyond the scope of secondary school studies there is evidence that such systems have been assessed. For example the Victorian External VCE Specialist Mathematics Written Examination 2 (2008) contained coupled solutions for the interaction between rabbits and foxes. The students were required to draw phase planes from these solutions. In this same paper another question required students' to interpret a direction field for certain first order differential equations.

Mathematics is an ever evolving discipline. Although many mathematical models exist whose solutions are quite complex, this need not prevent students from exploring such ideas. The phase plane method gives the student another tool in their ever expanding 'tool box' to explore an exciting area of mathematics.

Acknowledgement

Special thanks must go to John C. Polking, Professor Emeritus, Rice University, Houston who gave his whole hearted permission to use the pplane software used to plot the phase plane portraits. A free-to-use Java version of this software is located at http://math.rice.edu/~dfield/dfpp.html

References

Edelstein-Keshet, L. (2005). Mathematical models in biology. New York: Siam.

Jordan, D. & Smith, P. (2007). Nonlinear ordinary differential equations: An introduction for scientists and engineers (4th ed.). New York: Oxford.

Liu, B., Zhang, Y., & Chen, L. (2005). The dynamical behaviors of a Lotka-Volterra predator-prey model concerning integrated pest management. Nonlinear analysis: Real world applications, 6(2), 227-243.

Mahaffy, J. (2001). Lotka-Volterra models. Retrieved from http://www- rohan.sdsu.edu/~jmahaffy/courses/f01/math636/lectures/lotka/qualde2.html

Mkhatshwa, T. & Mummert, A. (2010). Modeling super-spreading events for infectious diseases: Case study SARS. Accessed from http://arxiv.org/PS_cache/ arxiv/pdf/1007/1007.0908v2.pdf

Murray, J. (2002). Mathematical biology 1: An introduction (3rd ed.). New York: Springer.

VCAA. (2008). Victorian External VCE Specialist Mathematics Written Examination 2: 2008. Retrieved from http://www.vcaa.vic.edu.au/vcaa/vce/ studies/mathematics/specialist/pastexams/2008/2008specmaths2

Doug Braiden

John Paul College, Daisy Hill, Queensland

<dbraiden@jpc.qld.edu.au>
COPYRIGHT 2011 The Australian Association of Mathematics Teachers, Inc.
No portion of this article can be reproduced without the express written permission from the copyright holder.