The effects of spatially-constrained treatment regions upon a model of wombat mange
Ivy J. Hindle, Lawrence K. Forbes, Stephen J. Walters, Scott Carver
TL;DR
This study models how treating wombats in specific areas can reduce mange mites, but success depends on initial disease state and treatment location.
Contribution
The paper introduces a spatially constrained treatment model for managing Sarcoptes scabiei in wombats using non-linear PDEs.
Findings
Targeted treatment can reduce mite numbers in treatment regions and increase wombat populations.
Treatment effectiveness depends on initial disease state and treatment location within the wombat population.
Understanding host-pathogen dynamics and spatial terrain is crucial for successful treatment.
Abstract
The use of therapeutic agents is a critical option to manage wildlife disease, but their implementation is usually spatially constrained. We seek to expand knowledge around the effectiveness of management of environmentally-transmitted Sarcoptes scabiei on a host population, by studying the effect of a spatially constrained treatment regime on disease dynamics in the bare-nosed wombat Vombatus ursinus. A host population of wombats is modelled using a system of non-linear partial differential equations, a spatially-varying treatment regime is applied to this population and the dynamics are studied over a period of several years. Treatment could result in mite decrease within the treatment region, extending to a lesser degree outside, with significant increases in wombat population. However, the benefits of targeted treatment regions within an environment are shown to be dependent on…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9- —University of Tasmania
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsDermatological diseases and infestations · Bird parasitology and diseases · Study of Mite Species
Introduction
The effective management of infectious pathogens in free-ranging animals has many complex challenges, such as the identification of diseased individuals or the removal or targeted treatment of diseased individuals (Wobeser 2002). An option for disease management is the use of therapeutic agents, such as anti-parasitic drugs (Martin et al. 2019; Skerratt 2003; Death et al. 2011); however, the delivery of therapeutic agents to free-ranging hosts is challenging in nature. Consequently, treatment programs in wildlife are often spatially restricted to confined areas of specified sub-populations (Martin et al. 2019). Understanding how effective spatially restricted treatment regimes actually are, for host populations, is therefore valuable for in situ disease management practice.
The parasitic mite Sarcoptes scabiei, causing the disease sarcoptic mange, is one of the most host-generalist of mammalian parasites (Astorga et al. 2018). This parasite has been documented to infect approximately one hundred and fifty mammalian species across the globe (Pence and Ueckermann 2002), and continues to have an expanding host species range (Tompkins et al. 2015; Martin et al. 2018b). The mite burrows into the host’s skin, leading to a range of pathological outcomes, including pruritus, alopecia, secondary bacterial infections and emaciation (Arlian and Morgan 2017; Astorga et al. 2018). It is both an animal welfare issue and occasional conservation issue, causing epizootics and host population declines (Pence and Ueckermann 2002). Known as scabies when infecting a human host, it is also among the thirty most prevalent human diseases (estimated point prevalence two hundred million people) (World Health Organisation 2020) and is considered a Neglected Tropical Disease by the World Health Organisation (Karimkhani et al. 2017).
Sarcoptic mange is the most important disease impacting the bare-nosed wombat (Vombatus ursinus) (Martin et al. 2018a), with an estimated mortality rate of 100% for untreated individuals (Martin et al. 2018a; Skerratt 2003). The parasite is believed to have been introduced to Australia by European colonists and their domestic animals (Skerratt et al. 1998; Fraser et al. 2016). Wombat populations that are exposed to mange experience a variety of population outcomes, including both epizootics with population decline and endemic disease with relatively stable populations (Driessen et al. 2021; Skerratt et al. 1998; Martin et al. 2018a). Transmission of S. scabiei is understood to be environmental via burrows (Skerratt et al. 1998), between which solitary wombats move every 4–10 days, sharing asynchronously. In situ attempts at disease management in wombats is relatively common at the individual level (Skerratt 2003), and there is increasing interest about pathogen management at population scales. For example, a recent population scale treatment attempt during an epizootic showed short-term pathogen management was achievable (Martin et al. 2019). In most instances, population scale S. scabiei management attempts in wombats will be spatially restricted, so understanding how this restriction impacts spatial disease dynamics can inform management interventions.
In this study, we seek to characterize the effects of a realistic treatment regime on a population of bare-nosed wombats suffering from endemic S. scabiei infection. This treatment will be spatially confined to a sub-region of the host range. We aim to expand upon the spatio-temporally varying model proposed by Hindle et al. (2022), by introducing treatment for disease and exploring the long-term dynamics introduced by a short-term, spatially restricted treatment plan. A common treatment regime for wombats focusses on the repeated delivery of treatment [(e.g. at least 1 ml of moxidectin per 10 kg of animal weight (Martin et al. 2019)] over a defined period (e.g. weekly for 12 weeks), which may be targeted at infected individuals or indiscriminately at a population scale. This treatment is applied topically to the wombat, usually via a burrow flap (Martin et al. 2019). We focus on indiscriminate treatment at a population scale [consistent with Martin et al. (2019)]. Two differing sub-regions are considered, one situated at the centre of the host population range, and the other at one corner, with the intention of discovering if one or the other results in a sustained reduction in mite population density and an improved outcome for the wombats.
The wombat–mite model and the governing partial differential equations that describe it are introduced in Sect. 2. The steady-state populations possible in this model, for uniform treatment, are briefly reviewed in Sect. 3 and follow the developments in Beeton et al. (2019). It is found that up to four steady states are possible, although one of these is non-physical since it would involve negative populations; the remaining three consist of a total-extinction state, a mite-free state and an endemic state in which infected wombats and mites co-exist. Driving the system toward the mite-free state is an obvious intention of this study. We briefly consider a linearization of the governing equations in Sect. 4, because the governing equations in Sect. 2 are nonlinear and so can only be solved approximately with the aid of computers. This linearized solution therefore provides an important check on the accuracy of the computationally obtained nonlinear results. Nevertheless, the linearized solution has severe restrictions, since it is only valid when populations remain close to their assumed steady-state values, and so numerical solutions are needed for the fully nonlinear model outlined in Sect. 2, and our method for generating such solutions is discussed in Sect. 5. Results of this work are presented in Sect. 6; we first compare the linearized solution with the nonlinear computational results, and then investigate the effects of treatment regimes. Some final remarks in Sect. 7 conclude the paper.
The wombat–mite model
A conceptual model illustrating the wombat–mite system and its dynamics is illustrated in Fig. 1, with rate parameters for the system described in Table 1. These parameters will be assumed to be those used for the results in Sect. 6, unless indicated in the text. Fig. 1A modified version of the model developed in Hindle et al. (2022). Here the transmission method is environmental (indirect) and the population of the mites is modelled explicitly, to represent this indirect transmission. In this diagram, S represents the sub-population of susceptible wombats, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_L$$\end{document} are asymptomatic carriers of the disease, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_H$$\end{document} are symptomatic infected wombats, F represents the proportion of fomites existing in the environment, and R are the treated wombats. The quantity r is the treatment rate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_R$$\end{document} is the relapse rate at which treated wombats re-enter the susceptible group, b is the birth rate of wombats, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} is the death rate without infection of wombats, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _L$$\end{document} is the death rate due to asymptomatic infection, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _H$$\end{document} is the death rate due to symptomatic infection, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} is the rate at which susceptible wombats are infected, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} is the rate at which asymptomatic wombats become symptomatic, f is the rate that infected wombats drop fomites into the environment and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F$$\end{document} is the rate that fomites within the environment dieTable 1Rate parameters for the PDE (see Eqs. (1)–(5) and Fig. 1) describing transmission of S. scabiei among wombatsParameterSymbolValuesSourceWombat spatial diffusion rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} 1/(1.5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 365 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 0.5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 0.5) Banks et al. (2002)Wombat birth rateb1/(3 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 365) Beeton et al. (2019), Martin et al. (2019) and Triggs (2009)Non-linear disease spread rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} 0.01 Tamura et al. (2021)Wombat death rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} 1/(15 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} 365) Triggs (2009)Disease progression rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} 1/30 Skerratt (2003) and Martin et al. (2019)Exposed wombat death rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _L$$\end{document} 1/60 Martin et al. (2019)Infected wombat death rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _H$$\end{document} 1/60 Martin et al. (2019)Mite drop ratef1/5 Beeton et al. (2019)Mite death rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F$$\end{document} 1/19 Arlian et al. (1984) and Browne et al. (2021)Treatment rater0.9/7 Martin et al. (2019)Treatment relapse rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_R$$\end{document} 1/7 Beeton et al. (2019)All rates are in days, and the wombat spatial diffusion rate is in km/day
We create a system of nonlinear equations in two spatial dimensions and time, expanding upon the work of Hindle et al. (2022). The governing partial differential equations (PDEs) describing this system are:
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}{} & {} \frac{\partial S}{\partial t} = \sigma \nabla ^2 S + bN(1-N) - \frac{\beta SF}{1+F} - (\mu +r) S + k_R R \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}{} & {} \frac{\partial I_L}{\partial t} = \sigma \nabla ^2 I_L + \frac{\beta SF}{1+F} - (\mu +\gamma +\mu _{L}+r )I_L \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}{} & {} \frac{\partial I_H}{\partial t} = \sigma \nabla ^2 I_H + \gamma I_L - (\mu +\mu _H +r)I_H \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}{} & {} \frac{\partial R}{\partial t} = \sigma \nabla ^2 R + \left( S+I_L+I_H \right) r - (\mu +k_R) R \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}{} & {} \frac{\partial F}{\partial t} = \left( I_L+I_H \right) f - \mu _F F. \end{aligned}$$\end{document}Here, the quantities r and f are rate constants, as indicated in Fig. 1. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nabla ^2$$\end{document} operator used in the above equations is the two dimensional Laplacian operator
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \nabla ^2 \equiv \frac{\partial ^2}{\partial x^2} + \frac{\partial ^2}{\partial y^2} \end{aligned}$$\end{document}and represents the movement of wombats throughout the system. This movement is thus modelled as a diffusion process, with wombats spreading throughout the region as time progresses, governed by the spreading parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} in Table 1.
The governing equations (1)–(5) describe the movement of the wombat populations within the spatial region as well as their change over time. The wombat population is split into four sub-populations, following conventional SEIR models of disease spread (Kermack and McKendrick 1927). Thus the total wombat population is N, and we postulate
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} N = S + I_L + I_H + R, \end{aligned}$$\end{document}in which S denotes the sub-population of wombats that are susceptible to, but not yet infected by S. scabiei. In addition, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_L$$\end{document} represents the wombats that are infected by the disease, but are carrying a sufficiently low population of mites that they would not be seen as being infected in a standard field study; in an SEIR model, our sub-population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_L$$\end{document} would thus correspond to the “exposed” group. The sub-population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_H$$\end{document} represents those wombats carrying a high mite load and so show symptoms of the disease, and R is the “recovered” group that have been treated and are therefore temporarily immune from immediate further infection by S. scabiei. Finally, the population index of the fomites is labelled F and represents the mites existing on a host body, surviving within the bedding chamber of wombat burrows (Martin et al. 2018a; Beeton et al. 2019). This population is modelled as a proxy for indirect environmental transmission of sarcoptic mange in wombat populations.
The transmission terms in Eqs. (1), (2) at which the susceptible wombats enter the lightly-infected class are modelled with the frequency-dependent rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta F / (1+F)$$\end{document} . This is a self-limiting infection term, such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F / (1+F) \rightarrow 1$$\end{document} as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F \rightarrow \infty $$\end{document} , and is consistent with current understanding of mange transmission among wombats (Beeton et al. 2019; Martin et al. 2019; Tamura et al. 2021). The function is intended to reflect the environmental–transmission nature of S. scabiei in wombat populations, where an individual’s probability of exposure is related more to contact with contaminated burrow environments, rather than a direct encounter with another infected animal. Consequently, exposure occurs largely independently of wombat population density.
To approximate the physical conditions within which the wombats live, we consider a barrier around our spatial region of interest. For simplicity, this will be assumed to be a rectangular domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-L< x < L$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-B< y < B$$\end{document} , and the borders of this region will be supposed to prevent any wombat movement through them, although movement along the barrier will be allowed. These physical conditions that are being modelled can be considered to be a reserve that has been fenced in, or an area with natural barriers surrounding it, such as oceans, rivers, hills or dense forest. This reflects the actual situation for the wombats living at Narawntapu National Park (NNP), as detailed by Martin et al. (2018a). On the idealized rectangular boundaries considered in the present model, these restrictions are represented mathematically by the boundary conditions
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial S}{\partial x} = \frac{\partial I_L}{\partial x} = \frac{\partial I_H}{\partial x} = \frac{\partial R}{\partial x} = 0 \quad \textrm{on} \quad x = \pm L \end{aligned}$$\end{document}and
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial S}{\partial y} = \frac{\partial I_L}{\partial y} = \frac{\partial I_H}{\partial y} = \frac{\partial R}{\partial y} = 0 \quad \textrm{on} \quad y = \pm B. \end{aligned}$$\end{document}In order to complete the mathematical model of this situation, it is necessary to impose appropriate initial conditions. In this paper, we will consider conditions under which the system starts at one of its possible steady-state configurations (see Sect. 3), but mites are added into some portion of the physical-space domain. Furthermore, we wish to allow the possibility that our treatment rate may be a function both of space and time, so that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r \equiv r(x,y,t)$$\end{document} . This must be taken into account in the numerical solution technique (see Sect. 5), and physically represents the fact that treatment may only be offered over some fraction of the total region available to the wombats, and that the treatment is only able to be provided over short times. This reflects the actual treatment regime attempted for the wombats living in Narawntapu National Park, documented in Martin et al. (2019), in which treatment stations (burrow flaps) were refilled weekly; this treatment ran for 12 weeks in total, and the therapeutic agent (moxidectin) was observed to drop exponentially in efficacy as each week passed, since the insecticide was degraded by ultraviolet light. We discuss the results obtained with this present model, for short bursts of treatment applied periodically to the system, in Sect. 6.2.
Some comments are required, too, concerning the reasoning behind the parameter values assumed in this current model, as outlined in Table 1. Firstly, the wombat spatial diffusion rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} is largely due to females relocating after successfully rearing a joey, and they move approximately half a kilometre away (Banks et al. 2002). The wombat birth rate b is one joey per female per 1.5 years, and the death rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} is based on a wombat lifespan of 15 years. The transmission rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} is consistent with that assumed by Beeton et al. (2019), and the disease progression rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} is based on a 30 day window before symptoms become clinically visible ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_L$$\end{document} moves to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_H$$\end{document} ). The death rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _H$$\end{document} associated with mange is 60 days when showing clinical signs, and thus infection to mortality is about 90 days. The mite drop rate f is consistent with Beeton et al. (2019). To estimate a constant value for the treatment rate r, it was assumed that within the treatment region, 90% of the theoretical maximum treatment was delivered each 7 days (Martin et al. 2019). The treatment relapse rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_R$$\end{document} is based on the treatment on a wombat lasting 7 days before the animal could return to the susceptible class and potentially be re-infected by S. scabiei (Beeton et al. 2019; Martin et al. 2019), and the environmental mite death rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F$$\end{document} is for optimal environmental conditions, amounting to 19 days based on laboratory studies (Arlian et al. 1984; Browne et al. 2021).
The steady-state populations
The governing system of equations (1)–(5) in Sect. 2 allows several equilibrium population states to exist, in situations where the treatment rate r is constant across the entire solution domain. These were discussed briefly by Beeton et al. (2019). For completeness, and in view of their importance in the present model, they will be discussed briefly in this section.
Steady-state populations occur when the sub-populations S, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_L$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_H$$\end{document} , R and the mite density F are independent of spatial variables x and y and time t. The treatment rate r, too, must be assumed to be constant. Then all the derivatives in (1)–(5) are set to zero, leaving a system of five nonlinear algebraic equations for the equilibrium populations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( S_{eq}, I_{L,eq}, I_{H,eq}, R_{eq}, F_{eq} \right) $$\end{document} .
It then follows at once from (3) that
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} I_{L,eq} = \frac{\theta _H}{\gamma } I_{H,eq} \end{aligned}$$\end{document}and (5) similarly yields
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} F_{eq} = \frac{f}{\mu _F} \left( \frac{\theta _H + \gamma }{\gamma } \right) I_{H,eq} \equiv \zeta _H I_{H,eq}. \end{aligned}$$\end{document}In these two expressions, and elsewhere in this paper, it is convenient to define the intermediate constants
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta _L= & {} \mu + \mu _L + r \nonumber \\ \theta _H= & {} \mu + \mu _H + r. \end{aligned}$$\end{document}We have also created the constant
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \zeta _H = \frac{f \left( \theta _H + \gamma \right) }{\mu _F \gamma }. \end{aligned}$$\end{document}The steady-state form of Eq. (2) is simply
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \beta S_{eq} F_{eq} = \left( \theta _L + \gamma \right) \left( 1 + F_{eq} \right) I_{L,eq} \end{aligned}$$\end{document}and when (9), (10) are used in this expression, we obtain
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \beta \gamma \zeta _H S_{eq} I_{H,eq} = \theta _H \left( \theta _L + \gamma \right) \left( 1 + \zeta _H I_{H,eq} \right) I_{H,eq}. \end{aligned}$$\end{document}This equation presents two possible solutions; either
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} I_{H,eq} = 0 \end{aligned}$$\end{document}or else
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} S_{eq} = \frac{\theta _H \left( \theta _L + \gamma \right) }{\beta \gamma } \left( \frac{1}{\zeta _H} + I_{H,eq} \right) . \end{aligned}$$\end{document}These two options are now studied using the remaining steady-state equations.
The simple choice in (13) is substituted into the steady-state forms of Eqs. (1) and (4), where they are easily seen to reveal two equilibrium solutions. The first of these is the total extinction case
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left( S_{eq}, I_{L,eq}, I_{H,eq}, R_{eq}, F_{eq} \right) = (0, 0, 0, 0, 0 ) \end{aligned}$$\end{document}with total wombat population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{eq} = 0$$\end{document} . The second equilibrium solution from this choice is seen to be
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}{} & {} I_{L,eq} = I_{H,eq} = F = 0; \nonumber \\{} & {} S_{eq} = \frac{b - \mu }{b K}; \quad R_{eq} = \frac{r(b - \mu )}{\left( \mu + k_R \right) b K}; \quad N_{eq} = \frac{b - \mu }{b}. \end{aligned}$$\end{document}Here, we have defined a further intermediate constant
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} K = \frac{\mu + k_R + r}{\mu + k_R} \end{aligned}$$\end{document}to simplify the notation. This second equilibrium (16) is surely the steady-state outcome to be desired, since it represents a situation in which the wombats live their lives free from mites.
The other choice (14) can give rise to circumstances in which both wombats and mites co-exist. It therefore needs to be explored, although the algebra is a little tedious. The steady-state (derivative-free) forms of the model equations (2), (3), (5) give
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}{} & {} S_{eq} = \mathcal {M} \left( 1 + \zeta _H I_{H,eq} \right) ; \quad { I_{L,eq} = \frac{\theta _H}{\gamma } I_{H,eq}; \quad } \nonumber \\{} & {} F_{eq} = \zeta _H I_{H,eq} {; \quad R_{eq} = \frac{r}{\mu + k_R} \left( \mathcal {M} + \mathcal {Z} I_{H,eq} \right) } \end{aligned}$$\end{document}and it follows from (6) that the total wombat number can be expressed as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} N = K \left( \mathcal {M} + \mathcal {Z} I_{H,eq} \right) . \end{aligned}$$\end{document}In Eqs. (18) and (19) we have defined further intermediate constants
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathcal {M} = \frac{\theta _H \left( \theta _L + \gamma \right) }{\beta \gamma \zeta _H}; \quad { \mathcal {Z} = \frac{\beta \left( \theta _H + \gamma \right) + \theta _H \left( \theta _L + \gamma \right) }{\beta \gamma } } \end{aligned}$$\end{document}and also made use of the constants defined in (11) and (12). These expressions (18) and (19) are substituted into the steady-state form of (1), and it is convenient to define a new quantity
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \xi = \mathcal {M} + \mathcal {Z} I_{H,eq} \end{aligned}$$\end{document}which may be shown to satisfy the equation
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} b K^2 \xi ^2 + \mathbb {T} \xi - \frac{\mathcal {M}}{\gamma \mathcal {Z}} \left( \theta _H \mu _L + \gamma \mu _H \right) = 0. \end{aligned}$$\end{document}Here, we have defined
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathbb {T} = -b K - \frac{k_R r}{\mu + k_R} + \frac{\zeta _H \mathcal {M}}{\mathcal {Z}} \left( \beta + \mu + r \right) \end{aligned}$$\end{document}and the constant K is as given in (17).
Since (22) is simply a quadratic equation, its two solutions can be written at once in the form
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\xi = \frac{1}{2b K^2} \left[ - \mathbb {T} \pm \sqrt{ \mathbb {T}^2 + 4b K^2 \frac{\mathcal {M}}{\gamma \mathcal {Z}} \left( \theta _H \mu _L + \gamma \mu _H \right) } \right] . } \end{aligned}$$\end{document}It is clear that both the solutions in (24) are real, and that the plus sign results in a positive value for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} and the negative sign always gives \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} negative. Once \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} has thus been determined, the equilibrium population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{H,eq}$$\end{document} is recovered from the expression (21), after which all the other populations can be obtained from (18) and (19).
The minus sign option in the solution (24) always gives \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi < 0$$\end{document} and therefore always results in negative values for the sub-population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{H,eq}$$\end{document} in (21). Consequently, this solution has no biological significance and can be discarded. However, the other solution, obtained from the choice of the plus sign in (24), is able to give positive and biologically relevant values for all the populations in this model, for appropriate choices of the parameters. This shows that a steady state is indeed possible, in which mites are endemic in the system. The aim of our treatment regime in such circumstances is to effect a movement from this endemic steady state to the mite-free state given in (16).Fig. 2. Equilibrium steady-state: a recovered sub-population R, b fomite population F and c total wombat population N, obtained using the parameters in Table 1. Results are shown for the mite death rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 1/19 \approx 0.0526$$\end{document} given in the Table, and also for the smaller value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 0.005$$\end{document} discussed in the text
Figure 2 shows how some steady-state populations are affected by the treatment rate r (here assumed constant over the entire domain). We show the recovered sub-population R of wombats, the population density F of mites and the total wombat population N in parts (a), (b) and (c) respectively. On the left of each figure, the mite death rate has been set to the value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 1/19 \approx 0.0526$$\end{document} given in Table 1, and the diagrams on the right have used the smaller value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 0.005$$\end{document} as a comparison. The solid line (blue online) represents the total extinction equilibrium (15) and the other (black) solid line corresponds to the mite-free equilibrium (16). The two dashed lines come from the solution (24), and these represent the two endemic states. The formula with the minus sign in (24) gives negative mite and wombat populations, as discussed above, and these have no biological meaning, although for completeness we have nevertheless drawn those unphysical negative populations in Fig. 2. We observe, too, that for the steady-state populations obtained with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 1/19 \approx 0.0526$$\end{document} on the left-hand side, the positive sign in the solution (24) gives an endemic equilibrium in which the recovered population R and total population N are both above the carrying capacity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=1$$\end{document} , when the treatment rate r is sufficiently large; however, such an unphysical situation would not occur, since in those cases the mite population F has fallen below zero. Instead, the solution would revert to the mite-free equilibrium (16) for these larger treatment rates, exactly as desired. In addition, it is interesting to note that this “plus” endemic equilibrium, under our parameter values and for the case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 1/19 \approx 0.0526$$\end{document} , is highly sensitive (unstable) to the introduction of even small amounts of treatment.
The smaller mite death-rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 0.005$$\end{document} on the right-hand side of each diagram in Fig. 2 is almost certainly unphysical, since it would represent free mites living in the environment for 200 days, and this has not been encountered in our observations. In spite of the fact that it is therefore not of biological interest, it has nevertheless been included here so as to illustrate the mathematical dynamics possible in a model such as this. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 0.005$$\end{document} with treatment rates in about the interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0< r < 0.1$$\end{document} , it is clear that one of three possible steady states could occur. The populations could fall to extinction (15), and the mite-free population (16) is at least a mathematically possible outcome; the stable state, however, is the endemic one indicated by the (red) dashed line, in which mites and wombats co-exist in this interval of treatment rates r. This corresponds to the “plus” sign in (24), which the linearized solution of the next Sect. 4 shows to be the stable one for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0< r < 0.1$$\end{document} . It is interesting to observe, too, that the mite population F actually first increases with increasing treatment rate r, as shown with the (red) dashed line in part (b) for mite death rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu =0.005$$\end{document} . This somewhat paradoxical situation comes about because, with low treatment rates r, the number of wombats increases slightly, making more hosts available for mite infestation. In particular, the numbers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_L$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_H$$\end{document} of exposed and infected wombats increase, and it follows from (5) that the growth rate of mites increases even as the wombat numbers start to recover. Eventually, however, as the treatment rate r is increased still further, the mite numbers do fall away below zero, at which point the mite-free equilibrium (16) becomes the state that is achieved.
A linearized approximation
In this section, we develop an approximation to the full system of Eqs. (1)–(5) presented in Sect. 2, based on an assumption that the populations never stray far from (one of) their equilibrium values calculated in Sect. 3. This will give rise to a system of linear partial differential equations, which will be solved essentially exactly, although computer evaluation of the final solution will still be required. The reason for creating and solving this approximate model is primarily to serve as an important check on the accuracy of the numerical solutions to the full nonlinear model (1)–(5) that are needed later in the paper. Thus the linearized solution of this section is a necessary component in the proper analysis of our nonlinear model; nevertheless, it is quite limited in its scope, because it is only able to assess the behaviour of the system close to one of its equilibria.
It will be supposed here that a steady-state situation has been achieved, in which the populations of the fomites, as well as each of the sub-populations of the wombats, has reached equilibrium, as in Sect. 3. This equilibrium state is then disturbed, and here, it will be assumed that the perturbation remains small in comparison to the equilibrium population sizes. The size of this disturbance is estimated by some dimensionless parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} that measures the change in a population from its equilibrium value. Each population S(x, y, t) and so on is then represented mathematically as a small perturbation to its equilibrium level, so that, overall
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} S(x,y,t)= & {} S_{eq} + \varepsilon S_1 (x,y,t) + \mathcal {O}(\varepsilon ^2) \nonumber \\ I_L (x,y,t)= & {} I_{L,eq} + \varepsilon I_{L1} (x,y,t) + \mathcal {O}(\varepsilon ^2) \nonumber \\ I_H (x,y,t)= & {} I_{H,eq} + \varepsilon I_{H1} (x,y,t) + \mathcal {O}(\varepsilon ^2) \nonumber \\ R(x,y,t)= & {} R_{eq} + \varepsilon R_1 (x,y,t) + \mathcal {O}(\varepsilon ^2) \nonumber \\ F(x,y,t)= & {} F_{eq} + \varepsilon F_1 (x,y,t) + \mathcal {O}(\varepsilon ^2). \end{aligned}$$\end{document}The task is now to find each of the perturbation functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_1 (x,y,t)$$\end{document} and so on in this expansion. If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_1$$\end{document} grows larger as time increases, then the solution S moves further away from its equilibrium population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{eq}$$\end{document} and as a result, that equilibrium is unstable. If, however, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_1$$\end{document} decreases with time, then the solution S approaches its equilibrium value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{eq}$$\end{document} which consequently is stable.
These perturbed expressions in the representation (25) are substituted into the fully nonlinear system (1)–(5) and only terms of the first power in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} are retained. This process is described by Hindle et al. (2022, Appendix A) and, after some algebra, yields
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial S_1}{\partial t}= & {} \sigma \nabla ^2 S_1 + b N_1 \left( 1 - 2 N_{eq} \right) - \frac{\beta F_{eq} S_1}{1 + F_{eq}} - \frac{\beta S_{eq} F_1}{\left( 1 + F_{eq} \right) ^2} -(\mu + r) S_1 + k_R R_1 \nonumber \\ \frac{\partial I_{L1}}{\partial t}= & {} \sigma \nabla ^2 I_{L1} + \frac{\beta F_{eq} S_1}{1 + F_{eq}} + \frac{\beta S_{eq} F_1}{\left( 1 + F_{eq} \right) ^2} - \left( \theta _L + \gamma \right) I_{L1} \nonumber \\ \frac{\partial I_{H1}}{\partial t}= & {} \sigma \nabla ^2 I_{H1} + \gamma I_{L1} - \theta _{H} I_{H1} \nonumber \\ \frac{\partial R_1}{\partial t}= & {} \sigma \nabla ^2 R_1 + \left( S_1 + I_{L1} + I_{H1} \right) r - \left( \mu + k_R \right) R_1 \nonumber \\ \frac{\partial F_1}{\partial t}= & {} \left( I_{L1} + I_{H1} \right) f - \mu _F F_1 \end{aligned}$$\end{document}in which we have made use of the two constants in (11).
Since these Eq. (26) are to be solved in the rectangular domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-L< x < L$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-B< y < B$$\end{document} , subject to Neumann boundary conditions of the form in (7), (8), their solutions are sought in the Fourier-series forms
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} S_1 (x,y,t)= & {} \sum _{m=0}^{\infty } \sum _{n=0}^{\infty } S^L_{m,n} (t) \cos \left( \frac{m\pi (x-L)}{2L} \right) \cos \left( \frac{n\pi (y-B)}{2B} \right) \nonumber \\ I_{L1} (x,y,t)= & {} \sum _{m=0}^{\infty } \sum _{n=0}^{\infty } E^L_{m,n} (t) \cos \left( \frac{m\pi (x-L)}{2L} \right) \cos \left( \frac{n\pi (y-B)}{2B} \right) \nonumber \\ I_{H1} (x,y,t)= & {} \sum _{m=0}^{\infty } \sum _{n=0}^{\infty } H^L_{m,n} (t) \cos \left( \frac{m\pi (x-L)}{2L} \right) \cos \left( \frac{n\pi (y-B)}{2B} \right) \nonumber \\ R_1 (x,y,t)= & {} \sum _{m=0}^{\infty } \sum _{n=0}^{\infty } R^L_{m,n} (t) \cos \left( \frac{m\pi (x-L)}{2L} \right) \cos \left( \frac{n\pi (y-B)}{2B} \right) \nonumber \\ F_1 (x,y,t)= & {} \sum _{m=0}^{\infty } \sum _{n=0}^{\infty } F^L_{m,n} (t) \cos \left( \frac{m\pi (x-L)}{2L} \right) \cos \left( \frac{n\pi (y-B)}{2B} \right) . \end{aligned}$$\end{document}The aim is now to calculate the time-dependent Fourier coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S^L_{m,n} (t)$$\end{document} , and so on, for this linearized approximation.
The assumed forms (27) for the solution are substituted into the linearized equations (26). Each equation is then multiplied by a typical basis function
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \cos \left( k\pi (x-L) / (2L) \right) \cos \left( \ell \pi (y-B) / (2B) \right) , \end{aligned}$$\end{document}for integers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k = 0, 1, 2, \dots $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell = 0, 1, 2, \dots $$\end{document} , and integrated over the rectangular solution domain, making use of the orthogonality of the trigonometric functions. Since the partial differential equations are linear and have constant coefficients, each Fourier mode \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(k, \ell )$$\end{document} de-couples from all the others, and after some algebra, the Fourier-analyzed partial differential equations (26) yield a system of five ordinary differential equations of the form
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\textrm{d} \textbf{V}^L_{k,\ell }}{\textrm{d} t} = \textbf{J}_{k,\ell } \textbf{V}^L_{k,\ell } \end{aligned}$$\end{document}independently for each mode \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(k,\ell )$$\end{document} . The system involves the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(5 \times 1)$$\end{document} vector
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{V}^L_{k,\ell } = \bigl [ S^L_{k,\ell } (t); E^L_{k,\ell } (t); H^L_{k,\ell } (t); R^L_{k,\ell } (t); F^L_{k,\ell } (t) \bigr ]^T \end{aligned}$$\end{document}which consists of all the time-dependent Fourier coefficients that have been defined in the assumed solution (27). It also involves the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(5 \times 5)$$\end{document} Jacobian matrix
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{J}_{k,\ell } = \left[ \begin{array}{c c c c c} -A_{k,\ell } &{}\quad \Gamma _{eq} &{}\quad \Gamma _{eq} &{}\quad \Gamma _{eq}+k_R &{}\quad -Q_{eq} \\ P_{eq} &{}\quad -B_{k,\ell } &{}\quad 0 &{}\quad 0 &{}\quad Q_{eq} \\ 0 &{}\quad \gamma &{}\quad -C_{k,\ell } &{}\quad 0 &{}\quad 0 \\ r &{}\quad r &{}\quad r &{}\quad -D_{k,\ell } &{}\quad 0 \\ 0 &{}\quad f &{}\quad f &{}\quad 0 &{}\quad -\mu _F \end{array} \right] \end{aligned}$$\end{document}in which, for ease of display, we have introduced seven additional intermediate constants. These are
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Gamma _{eq} = b \left( 1 - 2 N_{eq} \right) ;\quad P_{eq} = \frac{\beta F_{eq}}{1 + F_{eq}};\quad Q_{eq} = \frac{\beta S_{eq}}{\left( 1 + F_{eq} \right) ^2} \end{aligned}$$\end{document}and
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} A_{k,\ell }= & {} \sigma \Delta _{k,\ell }^2 - \Gamma _{eq} + P_{eq} + \mu + r \nonumber \\ B_{k,\ell }= & {} \sigma \Delta _{k,\ell }^2 + \theta _L + \gamma \nonumber \\ C_{k,\ell }= & {} \sigma \Delta _{k,\ell }^2 + \theta _H \nonumber \\ D_{k,\ell }= & {} \sigma \Delta _{k,\ell }^2 + \mu + k_R. \end{aligned}$$\end{document}The extra constant
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Delta _{k,\ell }^2 = \left( \frac{k\pi }{2L} \right) ^2 + \left( \frac{\ell \pi }{2B} \right) ^2 \end{aligned}$$\end{document}used in (33) is an effective squared total wavenumber for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(k,\ell )$$\end{document} Fourier mode.
The solution of the vector differential equation (29) can be written down at once. It is
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{V}_{k,\ell }^L (t) = \sum _{q=1}^5 \mathcal {C}_{k,\ell }^{(q)} \exp \left[ \lambda _{k,\ell }^{(q)} t \right] \textbf{v}_{k,\ell }^{(q)} \end{aligned}$$\end{document}in which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{k,\ell }^{(q)}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{v}_{k,\ell }^{(q)}$$\end{document} are the five eigenvalues and their corresponding eigenvectors, for the Jacobian matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{J}_{k,\ell }$$\end{document} in (31), at each of the Fourier modes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(k,\ell )$$\end{document} . The five constants \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {C}_{k,\ell }^{(q)}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 1, \dots , 5$$\end{document} , are so far arbitrary, but their values are determined from the initial values of the Fourier coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S^L_{k,\ell } (0)$$\end{document} and so on in the vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{V}^L_{k,\ell } (0)$$\end{document} in (30). If, for each Fourier mode \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(k,\ell )$$\end{document} , these initial values are given, then the five constants are found by solving the matrix equation
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{Q}_{k,\ell } \textbf{C}_{k,\ell } = \textbf{V}^L_{k,\ell } (0) \end{aligned}$$\end{document}for the vector
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{C}_{k,\ell } = \left[ \mathcal {C}_{k,\ell }^{(1)}; \dots ; \mathcal {C}_{k,\ell }^{(5)} \right] ^T. \end{aligned}$$\end{document}The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(5 \times 5)$$\end{document} matrix
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{Q}_{k,\ell } = \left[ \textbf{v}_{k,\ell }^{(1)}; \dots ; \textbf{v}_{k,\ell }^{(5)} \right] \end{aligned}$$\end{document}in the initial condition (36) is made up of the five eigenvectors of the Jacobian matrix (31).
In principle, this analysis now gives the full linearized solution in some neighbourhood of one of the equilibrium points in Sect. 3. All that is required to evaluate the solution is to specify the initial condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{V}^L_{k,\ell } (0)$$\end{document} in (30), for each Fourier mode. The eigenvalues \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{k,\ell }^{(q)}$$\end{document} then determine whether the equilibrium state is stable and small perturbations to equilibrium conditions die away with time, or unstable such that the linearized solution grows exponentially with time.
Calculating these eigenvalues for general parameter values is, however, extremely difficult. Consequently, it is often not possible to give precise conditions under which the three feasible steady states (15), (16), and (24) with the plus sign chosen, are stable or unstable. In those cases, the eigenvalues and eigenvectors must be evaluated numerically, and the solutions (35) then obtained for Fourier modes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k, \ell = 0, 1, 2, \dots $$\end{document} . Finally, the linearized population distributions are evaluated approximately from the Fourier series (27).
In spite of these difficulties, however, it turns out that some important analytical information is nevertheless able to be obtained about the stability of these steady states. Firstly, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(5 \times 5)$$\end{document} Jacobian matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{J}_{k,\ell }$$\end{document} in (31) simplifies very considerably when applied to the total extinction steady state (15), and its eigenvalues can, in fact, be calculated in closed form. After some algebra, they turn out to be
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}{} & {} \lambda _{k,\ell }^{(1)} = - \mu _F; \quad \lambda _{k,\ell }^{(2)} = - B_{k,\ell }; \quad \lambda _{k,\ell }^{(3)} = - C_{k,\ell }; \nonumber \\{} & {} \lambda _{k,\ell }^{(4)} = - \sigma \Delta _{k,\ell }^2 + \left( b - \mu \right) ; \quad \lambda _{k,\ell }^{(5)} = - \sigma \Delta _{k,\ell }^2 - \left( \mu + k_R + r \right) . \end{aligned}$$\end{document}The auxiliary constants in these expressions are as defined in (33) and (34). Now for high modes, in which integers k and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell $$\end{document} are arbitrarily large, these eigenvalues (37) are all real and negative, and so the high Fourier modes are all stable. The critical case is the zeroth mode \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k = \ell = 0$$\end{document} corresponding to a spatially uniform disturbance, when the fourth eigenvalue in the system (37) becomes
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda _{0,0}^{(4)} = \left( b - \mu \right) . \end{aligned}$$\end{document}This is positive when birth rate b exceeds the death rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} , so this zeroth mode is unstable. Consequently, the total extinction steady state (15) is unstable.
The mite-free steady state (16) is obviously the desired equilibrium configuration in this model, and so it is important to know under what conditions the linearized solution would predict it to be stable. Unfortunately, this information seems too difficult to obtain exactly, for general parameter values. However, some important information about its stability can be derived analytically in the special case when the additional death rates for the lightly and heavily infected wombats are equal, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _L = \mu _H$$\end{document} . In this case, it follows from (11) that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _L = \theta _H$$\end{document} . Then the eigenvalues of the Jacobian matrix (31) evaluated for the mite-free case (16) can, in fact, be calculated exactly. A similar situation was encountered by Beeton et al. (2019). After considerable algebra,1 the five eigenvalues of (31) for the mite-free state with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _L = \mu _H$$\end{document} are found to be
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda _{k,\ell }^{(1)}= & {} - \sigma \Delta _{k,\ell }^2 - \gamma - \theta _H \nonumber \\ \lambda _{k,\ell }^{(2)}= & {} - \sigma \Delta _{k,\ell }^2 - \left( b - \mu \right) \nonumber \\ \lambda _{k,\ell }^{(3)}= & {} - \sigma \Delta _{k,\ell }^2 - \mu - k_R - r \nonumber \\ \lambda _{k,\ell }^{(4,5)}= & {} \frac{1}{2} \left\{ - \left( C_{k,\ell } + \mu _F \right) \pm \sqrt{ \left( C_{k,\ell } + \mu _F \right) ^2 - 4 \left[ \mu _F C_{k,\ell } - \frac{f \beta (b - \mu )}{b K} \right] } \right\} . \nonumber \\ \end{aligned}$$\end{document}For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b > \mu $$\end{document} , the first three eigenvalues in (38) are all negative and the remaining two eigenvalues \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{k,\ell }^{(4,5)}$$\end{document} are both real. This mite-free equilibrium, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _L = \mu _H$$\end{document} is stable if and only if all these eigenvalues are negative, and this is guaranteed if
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mu _F C_{k,\ell } - \frac{f \beta (b - \mu )}{b K} > 0 \end{aligned}$$\end{document}for all spatial Fourier modes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k, \ell = 0, 1, 2, \dots $$\end{document} . This inequality is trivially true for high modes k and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell $$\end{document} , but again, the zeroth modes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k = \ell = 0$$\end{document} pose the greatest challenge. Setting both k and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell $$\end{document} to zero in this inequality and rearranging then gives the stability result
Theorem 1
If the death rates are equal across the two infected classes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_L$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_H$$\end{document} , that is, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _L = \mu _H$$\end{document} , the mite-free steady state (16) is stable if and only if
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mu _F > \frac{f \beta \left( b- \mu \right) \left( \mu + k_R \right) }{b \left( \mu + k_R + r \right) \left( \mu + \mu _H + r \right) }. \end{aligned}$$\end{document}This condition (39) was found also by Beeton et al. (2019) for the purely time-dependent model. Although this stability condition (39) is only strictly true for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _L = \mu _H$$\end{document} , we anticipate that it will give at least an approximate guide to the stability of this key mite-free state more generally, and so the region of stability predicted by (39) is shown in Fig. 3 for parameter values taken from Table 1. When stability condition (39) does not hold, then the endemic state (24) will become the stable situation.Fig. 3. Region of stability for the mite-free equilibrium, as predicted by the linearized solution for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _L = \mu _H$$\end{document} and parameter values taken from Table 1
To determine stability of all the equilibrium states more generally, including cases when Theorem 1 does not apply, the linearized solution and its eigenvalues must be determined numerically, and this is now briefly outlined. To fix ideas, let us suppose that all the populations are at their steady-state values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( S_{eq}, I_{L,eq}, I_{H,eq}, R_{eq}, N_{eq} \right) $$\end{document} . Then, at the initial time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0$$\end{document} , mites are suddenly introduced over some rectangular region within the overall domain; let us suppose that the centre of this mite-infested region is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(x,y) = \left( x_C, y_C \right) $$\end{document} and the region has size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2 x_R$$\end{document} by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2 y_R$$\end{document} . Therefore the region newly infested with mites lies over the rectangle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_C - x_R< x < x_C + x_R$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_C - y_R< y < y_C + y_R$$\end{document} and to ensure this lies wholly within our overall domain, we impose conditions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_R - L< x_C < L - x_R$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_R - B< y_C < B - y_R$$\end{document} . For simplicity, we will assume that the mite population within the infected rectangular region has the constant value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} . Consequently, the initial condition on the mites is
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} F(x,y,0) = F_{eq} + {\left\{ \begin{array}{ll} \varepsilon , &{} \quad \text {inside rectangle} \\ 0, &{} \quad \text {outside rectangle}. \end{array}\right. } \end{aligned}$$\end{document}(We observe that Eq. (40) now serves to define the small parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} .) This condition may now be used to determine the initial values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F^L_{m,n} (0)$$\end{document} of the Fourier series for the mite numbers, in the expression (27). This series representation is subject to Fourier analysis, by multiplying by the basis functions (28) and integrating over the entire domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-L< x < L$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-B< y < B$$\end{document} . A somewhat lengthy calculation yields
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} F^L_{0,0} (0)= & {} \left( \frac{x_R}{L} \right) \left( \frac{y_R}{B} \right) \nonumber \\ F^L_{0,\ell } (0)= & {} \frac{4}{\ell \pi } \left( \frac{x_R}{L} \right) \cos \left( \frac{\ell \pi \left( B - y_C \right) }{2B} \right) \sin \left( \frac{\ell \pi y_R}{2B} \right) , \quad \ell > 0 \end{aligned}$$\end{document}for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m=0$$\end{document} modes, along with
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} F^L_{k,0} (0) = \frac{4}{k\pi } \left( \frac{y_R}{B} \right) \cos \left( \frac{k\pi \left( L - x_C \right) }{2L} \right) \sin \left( \frac{k\pi x_R}{2L} \right) \quad k > 0 \end{aligned}$$\end{document}for the modes with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=0$$\end{document} . Finally, the general case gives
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} F^L_{k,\ell } (0)&{=}&\frac{16}{k\ell \pi ^2} \cos \left( \frac{k\pi \left( L {-} x_C \right) }{2L} \right) \sin \left( \frac{k\pi x_R}{2L} \right) \cos \left( \frac{\ell \pi \left( B {-} y_C \right) }{2B} \right) \sin \left( \frac{\ell \pi y_R}{2B} \right) \nonumber \\{} & {} \quad \text {for } \quad k = 1, 2, 3, \dots , \quad \ell = 1, 2, 3, \dots . \end{aligned}$$\end{document}Because the solution is starting from an equilibrium population that has been perturbed by the introduction of mites at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0$$\end{document} , the initial vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{V}^L_{k,\ell }$$\end{document} in (30) becomes simply
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{V}^L_{k,\ell } = \bigl [ 0; 0; 0; 0; F^L_{k,\ell } (0) \bigr ]^T. \end{aligned}$$\end{document}Furthermore, because the approximate system in this Section is linear (with constant coefficients), then each Fourier mode \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( k, \ell \right) $$\end{document} acts completely independently of all the other modes, as discussed above. Consequently, at each separate mode \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( k, \ell \right) $$\end{document} , the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$( 5 \times 1)$$\end{document} matrix equation (36) can be solved for the five constants in the vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{C}_{k,\ell }$$\end{document} using the initial conditions (44), with values for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F^L_{k,\ell } (0)$$\end{document} taken from whichever of the expressions (41)–(43) is appropriate for the Fourier mode under consideration. The solution (27) can now be evaluated at any desired time, using the formula (35) for its time-dependent Fourier coefficients.
In practice, we choose to smooth the discontinuous initial condition (40), since it is well known that Fourier series of the type (27) cannot converge near a discontinuity; instead, the reconstructed function contains undesirable oscillations in that vicinity. This is known as the Gibbs phenomenon (Kreyszig 2011, page 515), and can be a particular problem in population models in which a population drops to zero, since spurious small oscillations might then produce negative population values instead of exactly zero, and this of course is unphysical. To avoid this, we use Lanczos smoothing (see Duchon 1979). To illustrate this for a single independent variable x, consider a periodic function f(x) over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-L< x < L$$\end{document} represented spectrally as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(x) = \sum _{m=0}^{\infty } A_m \cos \left( m\pi (x-L) / (2L) \right) . \end{aligned}$$\end{document}If f(x) has a discontinuity at some point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_D \in (-L, L)$$\end{document} , this series will produce oscillations near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_D$$\end{document} according to the Gibbs phenomenon. To avoid this, the Fourier coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_m$$\end{document} are replaced with new coefficients
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \overline{A}_m = A_m \frac{\sin \left( \sigma _L m \right) }{\sigma _L m}. \end{aligned}$$\end{document}This replaces the original discontinuous function f(x) with a new function which can be shown to be the average of f(x) over a moving window \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( x-\sigma _L, x+\sigma _L \right) $$\end{document} centred at the point x. For appropriate choices of the Lanczos parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _L$$\end{document} , the new function so produced closely approximates the original function except that the discontinuity is replaced by a rapidly-varying smooth change without overshoot. In practice, we choose the Lanczos parameter to be about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _L \approx 0.05$$\end{document} .
Numerical solution of nonlinear problem
We now seek a numerical solution to the full system (1)–(5) of nonlinear partial differential equations in Sect. 2, that can describe large-amplitude deviations from the equilibrium populations discussed in Sect. 3. We make use of a spectral method, based on insights gained from the linearized solution (25), (27) in Sect. 4, and now seek a numerical solution to the nonlinear equations in the truncated Fourier-series form
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} S(x,y,t)= & {} \sum _{m=0}^M \sum _{n=0}^N S_{m,n} (t) \cos \left( \frac{m\pi (x-L)}{2L} \right) \cos \left( \frac{n\pi (y-B)}{2B} \right) \nonumber \\ I_L (x,y,t)= & {} \sum _{m=0}^M \sum _{n=0}^N E_{m,n} (t) \cos \left( \frac{m\pi (x-L)}{2L} \right) \cos \left( \frac{n\pi (y-B)}{2B} \right) \nonumber \\ I_H (x,y,t)= & {} \sum _{m=0}^M \sum _{n=0}^N H_{m,n} (t) \cos \left( \frac{m\pi (x-L)}{2L} \right) \cos \left( \frac{n\pi (y-B)}{2B} \right) \nonumber \\ R(x,y,t)= & {} \sum _{m=0}^M \sum _{n=0}^N R_{m,n} (t) \cos \left( \frac{m\pi (x-L)}{2L} \right) \cos \left( \frac{n\pi (y-B)}{2B} \right) \nonumber \\ F(x,y,t)= & {} \sum _{m=0}^M \sum _{n=0}^N F_{m,n} (t) \cos \left( \frac{m\pi (x-L)}{2L} \right) \cos \left( \frac{n\pi (y-B)}{2B} \right) \end{aligned}$$\end{document}In this representation, the integers M and N should be taken as large as possible. These forms are chosen to satisfy the boundary conditions (7) and (8) identically.
The aim of the numerical solution approach is to find the sets of Fourier coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{m,n} (t), \dots , F_{m,n} (t)$$\end{document} in the representation (45). To do this, these Fourier series are substituted into the full nonlinear set (1)–(5), then multiplied by the same sets of basis functions (28) that were used for the linearized solution in Sect. 4. In this case, however, the integers in these basis functions only take values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k = 0, 1, 2, \dots , M$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell = 0, 1, 2, \dots , N$$\end{document} , since the series in the representations (45) are truncated to total orders M and N. The resulting equations are now integrated over the solution domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-L<x < L$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-B< y < B$$\end{document} as before, making use of the orthogonality of the cosine functions.
Unlike the linearized solution of Sect. 4, where all the Fourier modes de-coupled and so could be treated independently of each other, the nonlinearity of the full system of equations causes each Fourier coefficient to be affected by all the others in the system. In addition, here we also want to allow the treatment rate r to be a function r(x, y, t) of both position and time. Consequently, the Fourier decomposition of the system (1)–(5) now gives rise to a large coupled nonlinear system, which we present here for completeness. The susceptible sub-population yields
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\textrm{d} S_{k,\ell } (t)}{\textrm{d} t}= & {} \frac{1}{\delta _{k,\ell } LB} \int _{-B}^{B} \int _{-L}^{L} \biggl [ bN(1 - N) - \frac{\beta SF}{1+F} - rS \biggr ] \nonumber \\{} & {} \quad \times \cos \left( \frac{k\pi (x-L)}{2L} \right) \cos \left( \frac{\ell \pi (y-B)}{2B} \right) \, \textrm{d}x \, \textrm{d}y \nonumber \\{} & {} \quad - \sigma \Delta _{k,\ell }^2 S_{k,\ell } (t) - \mu S_{k,\ell } (t) + k_R R_{k,\ell } (t) \end{aligned}$$\end{document}and the lightly infected (“exposed”) group results in the system of ordinary differential equations
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\textrm{d} E_{k,\ell } (t)}{\textrm{d} t}= & {} \frac{1}{\delta _{k,\ell } LB} \int _{-B}^{B} \int _{-L}^{L} \biggl [ \frac{\beta SF}{1+F} - r I_L \biggr ] \nonumber \\{} & {} \quad \times \cos \left( \frac{k\pi (x-L)}{2L} \right) \cos \left( \frac{\ell \pi (y-B)}{2B} \right) \, \textrm{d}x \, \textrm{d}y \nonumber \\{} & {} \quad - \sigma \Delta _{k,\ell }^2 E_{k,\ell } (t) - \left( \mu + \gamma + \mu _L \right) E_{k,\ell } (t). \end{aligned}$$\end{document}The heavily infected sub-population of wombats is modelled by (3), which after Fourier analysis gives rise to
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\textrm{d} H_{k,\ell } (t)}{\textrm{d} t}= & {} - \frac{1}{\delta _{k,\ell } LB} \int _{-B}^{B} \int _{-L}^{L} r I_H \cos \left( \frac{k\pi (x-L)}{2L} \right) \cos \left( \frac{\ell \pi (y-B)}{2B} \right) \, \textrm{d}x \, \textrm{d}y \\{} & {} \quad - \sigma \Delta _{k,\ell }^2 H_{k,\ell } (t) + \gamma E_{k,\ell } (t) - \left( \mu + \mu _H \right) H_{k,\ell } (t) \end{aligned}$$\end{document}and the equation for the recovered wombats yields a further system of differential equations
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\textrm{d} R_{k,\ell } (t)}{\textrm{d} t}= & {} \frac{1}{\delta _{k,\ell } LB} \int _{-B}^{B} \int _{-L}^{L} r \left[ S + I_L + I_H \right] \nonumber \\{} & {} \quad \times \cos \left( \frac{k\pi (x-L)}{2L} \right) \cos \left( \frac{\ell \pi (y-B)}{2B} \right) \, \textrm{d}x \, \textrm{d}y \nonumber \\{} & {} \quad - \sigma \Delta _{k,\ell }^2 R_{k,\ell } (t) - \left( \mu + k_R \right) R_{k,\ell } (t). \end{aligned}$$\end{document}The final partial differential equation (5) for the mite population results in the system
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\textrm{d} F_{k,\ell } (t)}{\textrm{d} t} = f \left[ E_{k,\ell } (t) + H_{k,\ell } (t) \right] - \mu _F F_{k,\ell } (t). \end{aligned}$$\end{document}These equations apply to all the Fourier modes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k = 0, 1, 2, \dots , M$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell = 0, 1, 2, \dots , N$$\end{document} in the truncated series representation (45).
The constants \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta _{k,\ell }^2$$\end{document} appearing in the system (46)–(49) are the same as those in (34) used in the linearized solution of Sect. 4. For convenience of notation, we have also defined constants
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \delta _{k,\ell } = \left\{ \begin{array}{ll} 4, &{}\quad \text {if } \quad k=0 \text { and } \ell = 0 \\ 2, &{}\quad \text {if either } \quad k=0 \text { or } \ell = 0 \\ 1, &{}\quad \text {if }\quad k \ne 0 \text { and } \ell \ne 0 \end{array} \right. \end{aligned}$$\end{document}in this governing system of equations for the Fourier coefficients.
These relations (46)–(49) constitute a system of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5(M+1)(N+1)$$\end{document} coupled nonlinear ordinary differential equations for the sets of coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{k,\ell }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{k,\ell }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_{k,\ell }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{k,\ell }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{k,\ell }$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k = 0, 1, 2, \dots , M$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell = 0, 1, 2, \dots , N$$\end{document} . We use the numerical integration routine ode45 provided in the programming language matlab to integrate these equations forward in time, starting from appropriate initial conditions, which are chosen to mimic those used for the linearized solution in Sect. 4. That is, all the coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{k,\ell } (0)$$\end{document} and so on are set to zero at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0$$\end{document} , so that the solution starts from a chosen equilibrium state, except for the coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{k,\ell } (0)$$\end{document} for the initial number density of mites. These are chosen to be the same as for the linearized solution, and are given by (41)–(43), so as to generate the initial rectangular deposition of mites described by (40), although in practice we again smooth this discontinuous initial mite distribution using the Lanczos filter, as discussed at the end of Sect. 4. Once these coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{k,\ell } (t)$$\end{document} and so on have been determined, by integration of the system (46)–(49), the populations of wombats and mites are then reconstructed from their Fourier-series representations (45).
Presentation of results
Comparison of linearized and nonlinear solutions
Fig. 4A comparison of results obtained from the linearized solution of Sect. 4 with those obtained using the algorithm in Sect. 5 for the nonlinear problem. The parameters are all as given in Table 1. This figure illustrates the mite numbers F at the initial time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0$$\end{document} and five subsequent times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 10$$\end{document} , 20, 30, 40 and 50, on the centre-plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y=0$$\end{document} . Linearized results are shown with solid (blue online) curves and numerical nonlinear results are drawn with dashed (red) lines (color figure online)
We begin this presentation of results by checking the agreement between the numerical algorithm of Sect. 5 with the exact solution for the linearized problem in Sect. 4. To do this, we take the parameters in the model to be those given in Table 1; the solution is started at the mite-free equilibrium (16), except that mites are added according to the initial condition (40) with maximum mite numbers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon = 0.2$$\end{document} . For simplicity, the rectangle over which the mites are initially introduced into the system is centred at the origin, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( x_C, y_C \right) = (0,0)$$\end{document} , with half-lengths \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_R = 1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_R = 1$$\end{document} kilometres. (Lanczos smoothing is then applied to the initial mites distribution, as described at the end of Sect. 4). This is a case in which good agreement between the linearized solution (Sect. 4) and the full nonlinear results (Sect. 5) is to be expected, as the mite population diffuses and dies, and the system returns to its mite-free state (16).
Figure 4 shows a centre-line profile of the mite population as a function of distance x, for the initial time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0$$\end{document} and a further five times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 10, 20, \dots , 50$$\end{document} days. These profiles were taken from the linearized solution (27) and the numerical nonlinear solution (45) by evaluating those expressions on the plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y = 0$$\end{document} . At initial time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0$$\end{document} , the centre-line profile for the mites consists simply of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon = 0.2$$\end{document} mites over the interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< x < 1$$\end{document} , but no mites otherwise, according to the initial condition (40), and this appears in the first diagram at the left of the figure. The nonlinear results presented in this diagram have all been obtained using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(M,N) = (101,101)$$\end{document} Fourier coefficients in the truncated series (45), and five times that number of numerical mesh points were used in the physical (x, y) space. This is done so as to ensure that the numerical evaluations of the integrals in the formulae (46)–(48) are all performed to sufficient accuracy. In Fig. 4, the linearized solution is shown with solid lines (blue online) and the nonlinear results are drawn using (red) dashed lines.
The agreement between the linearized and nonlinear solutions for this case is excellent. For the five sample times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 10, \dots , 50$$\end{document} illustrated in Fig. 4, there is no discernible difference between the two sets of results visible in the figure. This agreement continues to hold at later times than shown here, and has been checked carefully. This is an important step in the verification of the nonlinear solution algorithm in Sect. 5, and it gives confidence in our spectral method for solving the full nonlinear equations (1)–(5). For this example, the total wombat number started at the equilibrium value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{eq} = 0.8$$\end{document} in (16) and stayed very close to that value throughout, and so has not been presented here.Fig. 5A comparison of results obtained from the linearized solution of Sect. 4 with those obtained using the algorithm in Sect. 5 for the nonlinear problem. The parameters are as given in Table 1, except that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 0.005$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 0.05$$\end{document} . This figure illustrates: a the mite numbers F, and b the total wombat numbers N at the initial time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0$$\end{document} and five subsequent times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 200$$\end{document} , 400, 600, 800 and 1000, on the centre-plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y=0$$\end{document} . The linearized solution is drawn with solid (blue) lines and the nonlinear results with dashed (red) lines (color figure online)
In Fig. 5, we continue this comparison of the predictions of the approximate linearized theory (from Sect. 4) with the numerical solution of the fully nonlinear equations (using the scheme in Sect. 5). In this example, however, we now choose parameter values for which the linearized solution indicates that the “plus” steady state should be the stable one, and accordingly we choose this “plus” state (24) as the initial condition, except that the system is again also perturbed by adding a patch of mites over the same rectangular region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< x < 1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< y < 1$$\end{document} centred at the origin, as for Fig. 4. Again we assume initial disturbance magnitude \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon = 0.2$$\end{document} for the mite numbers in the rectangle.
For the example presented in Fig. 5, the mite death rate has been set at the (admittedly unrealistic) smaller value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 0.005$$\end{document} , so as to explore the dynamical behaviour possible in this system. The treatment rate is also reduced significantly to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 0.05$$\end{document} . It can be seen from the diagrams on the right-hand side of Fig. 2 that, at this treatment rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 0.05$$\end{document} , it is to be expected that the endemic equilibrium (24)—with the “plus” sign chosen for the square-root term—will now become the stable steady state, rather than the mite-free case (16). This is indeed confirmed by the linearized solution of Sect. 4 and the stability diagram Fig. 3 . This would be a less desirable situation from the biological point of view, even if it were achievable, since now both the wombats and the mites would survive in the environment, under these changed conditions.
The mite population F on the centreline \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y = 0$$\end{document} is shown in Fig. 5a and the total wombat numbers N are displayed in part (b). In each figure, the first diagram on the left gives the situation at the initial time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0$$\end{document} . Thus the total wombat numbers in part (b) are calculated from the endemic equilibrium (24) for the “plus” branch, giving \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{eq} = 0.3121$$\end{document} . The steady-state mite numbers with treatment rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 0.05$$\end{document} are calculated to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{eq} = 0.3337$$\end{document} and this likewise serves as the initial mite number density, except over the interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< x < 1$$\end{document} where additional mite numbers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon = 0.2$$\end{document} have been added, to give \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 0.5337$$\end{document} in that interval. This perturbed initial situation is shown in Fig. 5a. As expected, the perturbation to the mite numbers decays away as time proceeds, although the decay rate is clearly much slower than that shown in Fig. 4, and a much longer time interval is required, in order to see this effect. Thus, Fig. 5 gives mite and wombat numbers at the five later times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 200, 400, \dots , 1000$$\end{document} . Eventually, the mite numbers return to their steady-state value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{eq} = 0.3337$$\end{document} across the entire spatial domain and the wombat sub-populations also return slowly to their equilibrium values. It will be seen from Fig. 5 that, while there is generally good agreement between the linearized solution (sketched with solid lines) and the nonlinear solution (drawn with dashed lines), the agreement between them is not as precise as was the case in Fig. 4. This is not a numerical error on the part of the nonlinear solution, since convergence of the Fourier series in the representation (45) has been monitored very carefully; the results in Fig. 5 have been generated with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(M,N) = (61,61), (81,81)$$\end{document} and (101, 101) modes, and are indistinguishable from each other. Rather, the slight discrepancies between the linearized approximation and the nonlinear results, near the maxima or minima of the curves shown in Fig. 5, are an indication of the effects of nonlinearity on the spatial profiles, near these crests or troughs.Fig. 6A comparison of results obtained from the linearized solution of Sect. 4 with those obtained using the algorithm in Sect. 5 for the nonlinear problem. The parameters are as given in Table 1, except that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 0.005$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 0.05$$\end{document} . This figure illustrates: a the mite numbers F, and b the total wombat numbers N at times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=120$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 480$$\end{document} , 960, 1440, 1920 and 2400 days, on the centre-plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y=0$$\end{document} . The initial condition consisted of a perturbation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon =0.2$$\end{document} to the mite-free state (16). The linearized solution (dashed, red lines) soon fails, but the nonlinear solution (solid, blue lines) evidently makes a transition to the endemic “plus” steady state (color figure online)
So far, the linearized solution from Sect. 4 has given a reasonably good indication of the behaviour of the fully nonlinear wombat–mite system. This is because we have only used initial conditions that consist of a moderate perturbation to a stable equilibrium state; as time evolves, the perturbations die away and the system returns to equilibrium. This has given us an important check on the reliability of the nonlinear solution algorithm in Sect. 5.
In Fig. 6, however, this situation now changes radically. The parameter values used in this case are exactly the same as for Fig. 5; but the initial condition has been changed to become a perturbation of the mite-free equilibrium (16), rather than of the (stable) equilibrium (24), for the “plus” case. Consequently, the linearized solution now fails after a relatively short time. This is because linearization is only valid when all the populations remain close to an equilibrium state, as indicated in the expression (25). Linearized solutions are very limited in the behaviours they can exhibit, since, as (35) indicates, they can only either decay exponentially in time (to a stable equilibrium) or else grow exponentially (if the nearby equilibrium is unstable). In Fig. 6, the initial condition consisted of a perturbation to the mite-free equilibrium, which is unstable for these parameter values; as a result, the only behaviour the linearized solution can develop is for it to increase without bound as time evolves. This is evident in Fig. 6, where there is moderate agreement with the nonlinear solution at the earliest time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 120$$\end{document} shown, but by time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 480$$\end{document} the linearized solution is predicting an unrealistic explosion in mite numbers. We do not show the linearized solution for later times, since it ceases to be of any further value.
By contrast, the nonlinear solution in Fig. 6 reveals some interesting behaviour, of a fundamentally nonlinear nature. At the initial time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0$$\end{document} , the mite population is zero throughout the entire region, with the exception of the rectangle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1<x < 1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< y < 1$$\end{document} , where it has the perturbed value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon = 0.2$$\end{document} . By the earliest time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 120$$\end{document} shown in Fig. 6, the number of mites in this inner rectangle has grown slightly, as is evident in Fig. 6a. The (nonlinear) solution then grows much larger, and by day number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 960$$\end{document} in part (a), it has reached an approximate maximum of about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 1.5$$\end{document} at the centre of the region. After that time, the peak mite numbers begin to decrease; furthermore, the spatial pattern begins to diminish as the mite population distributes itself more and more uniformly across the environment, here represented by the rectangle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-L< x < L$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-B< y < B$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L = B = 5$$\end{document} kilometres. The mite population in Fig. 6 is evidently making a transition from the initial perturbed mite-free state (16) to the stable endemic “plus” equilibrium state (24), although it makes a very large amplitude excursion as it does so. The effect on the total wombat numbers in Fig. 6b is consistent with that interpretation of these results, since initially those numbers started at the mite-free value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{eq} = 0.8$$\end{document} given in (16), but by the last time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 2400$$\end{document} shown, they were evidently approaching the “plus” equilibrium value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{eq} = 0.3121$$\end{document} calculated from (24).Fig. 7. Some spatial patterns for selected times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=120$$\end{document} , 960 and 2400, for the same case illustrated in Fig. 6. This figure illustrates: a the mite numbers F, and b the total wombat numbers N
To illustrate further the role that spatial variability can have in the distribution of mites and wombats, we present in Fig. 7 some spatial patterns at the three selected times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 120$$\end{document} , 960 and 2400, for the same case as considered in Fig. 6. At the earliest time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 120$$\end{document} , the situation is still moderately close to the initial condition obtained from the mite-free equilibrium (16); thus the initial mite numbers are effectively still zero over the environment, except roughly within the rectangle into which they were introduced with mite density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = \varepsilon = 0.2$$\end{document} . The total wombat population N in part (b) is still close to the initial mite-free value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{eq} = 0.8$$\end{document} , although there is a small decrease in wombats within the rectangular region into which the mites were introduced.
The second time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 960$$\end{document} in Fig. 7 corresponds roughly to the moment at which mite numbers reach their maximum at the centre of the region, as shown in Fig. 6. Nevertheless, although mite numbers have risen to about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 1.5$$\end{document} at the centre of the region, at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 960$$\end{document} , they are nevertheless still almost zero at the corners \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(x,y) = ( \pm 5, \pm 5)$$\end{document} , as is evident at that time in Fig. 7a. The total wombat numbers in the centre of the region, shown in part (b), have decreased from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N = 0.8$$\end{document} to about 0.6. After this time the mite numbers decrease slowly, but the wombat numbers become more uniform in space, and decrease very significantly to the point that, by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 2400$$\end{document} days, they are approaching the endemic steady-state value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{eq} = 0.3121$$\end{document} .
Variable treatment rate—nonlinear results
For the example run shown in Figs. 6 and 7, the choice of parameters in the model ensured that the endemic “plus” equilibrium state (24) was the stable one. This means that, as time progresses, the populations within the environment eventually converge to that state. The nonlinear numerical method of Sect. 5 is capable of tracking those changes through time and space, even when the populations make a transition from one equilibrium state to another, as in the example shown. However, as with many dynamical systems of this type, the populations here may experience transitory large-amplitude excursions in which the populations are very far from any of the equilibrium states in Sect. 3, as they make the transition from their initial configuration to their eventual state. This was seen at about time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 960$$\end{document} in Figs. 6 and 7, for example. It is easy to imagine that, in the field, a wildlife manager may notice the development of spatial regions in which mite concentration becomes large, as in Fig. 6a, and may decide to intervene. For wombats in Narawntapu National Park (NNP) in Tasmania, this intervention is accomplished using an insecticide that is introduced onto the wombat skin by means of flaps at the entrances of their underground burrows, as documented by Martin et al. (2019). This treatment runs for about 12 weeks. It is probably not possible, in most cases, to treat the entire environment in this manner, and so the wildlife manager would quite reasonably focus resources on that portion where mite numbers are at their greatest. In this Section, we therefore seek to model how such a strategy alters mite and wombat numbers.
To implement an intervention strategy of this type within our mathematical model, it is now necessary to allow the treatment rate r to become a (known) function of both space and time. We do this here by postulating a spatio-temporal function of the type
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} r(x,y,t) = r_A + \left( r_T - r_A \right) f_T (t) g_T (x,y). \end{aligned}$$\end{document}Here, the constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A$$\end{document} represents a constant “ambient” treatment rate that is applied to the whole environment, perhaps by spraying with a mite growth-retarding hormone, for example. The additional function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_T (t) g_T (x,y)$$\end{document} is the extra treatment applied at specific times to targeted sub-regions in space. Since this is only applied over a sub-region within the overall domain, we model the spatial component of the treatment function as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} g_T (x,y) = {\left\{ \begin{array}{ll} 1, &{} \quad \text {inside rectangle } \mathcal {R}_T \\ 0, &{} \quad \text {outside rectangle } \mathcal {R}_T \end{array}\right. } \end{aligned}$$\end{document}similarly to the form (40) used for the initial condition. Here, the symbol \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_T$$\end{document} denotes some treatment region, and for simplicity we choose it to be the rectangular region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| x - x_0 | < x_T$$\end{document} ; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| y - y_0 | < y_T$$\end{document} with centre \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( x_0, y_0 \right) $$\end{document} and area \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4 x_T y_T $$\end{document} .
The choice \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 0.005$$\end{document} for the mite death rate used in Figs. 6 and 7 is not biologically reasonable, as observed above, since it would imply that free mites could survive for 200 days, and that is far longer than actually observed. This choice was made only to illustrate the dynamics that can be achieved by a model of this type, since it is possible that related models in slightly altered circumstances might also encounter such behaviour, but in biologically achievable situations [for example, the fungal pathogen causing Bat White Nose syndrome is long lived in the environment (Hoyt et al. 2021)]. For the remainder of this paper, we now focus on parameter values specifically relevant to wombat–mite interactions. These are summarized in Table 1.Fig. 8. Comparison of treated case (solid line—blue online) with untreated case (chain-dot line—purple online). Results obtained from the algorithm in Sect. 5 for the nonlinear problem. The parameters are as given in Table 1, except that there is assumed to be no ambient treatment, so that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A = 0$$\end{document} . There is then a single treatment (at the rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 0.9 / 7$$\end{document} given in Table 1), over the time interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$200< t < 250$$\end{document} , in the square \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< x,y < 1$$\end{document} centred at the origin. This figure illustrates: a the mite numbers F, and b the total wombat numbers N, at times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=120$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 240$$\end{document} , 720, 1320, 1920 and 2400 days, on the centre-plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y=0$$\end{document} . The initial condition consisted of a perturbation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon =0.2$$\end{document} to the mite-free state (16) (color figure online)Fig. 9A graph of the switching function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_T (t)$$\end{document} described in the text, that allows the treatment to be applied occasionally. Here, treatment is switched on at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=200$$\end{document} days and off again at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=250$$\end{document} days. This treatment regime is repeated every 250 days
Figure 8 shows a case in which the system is started from the mite-free equilibrium (16), but with an initial injection of mites over the unit square at the centre of the region, of perturbation magnitude \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon =0.2$$\end{document} . To ensure the accuracy of the results, the numbers of Fourier modes used were \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(M,N) = (121,121)$$\end{document} , with five times those numbers of mesh points. Here, however, a single treatment has also been applied at the centre of the physical region where the mite concentration is highest, perhaps mimicking the behaviour of a wildlife manager in the region. The treatment here has been applied over the 50-day interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$200< t < 250$$\end{document} , and over the square region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< x < 1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< y < 1$$\end{document} in the centre of the wombat environment. Thus the treatment zone \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_T$$\end{document} in (52) has parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( x_0, y_0 \right) = (0,0)$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_T = y_T = 1$$\end{document} . The model parameters are as given in Table 1. In addition, it is assumed in Fig. 8 that there is no ambient treatment over the general environment, so that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A = 0$$\end{document} , and the rate within the treatment region has been chosen to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_T = 0.9/7 = 0.1286$$\end{document} , as in Table 1. The solid (blue) lines in Fig. 8a, b show, respectively, the number F of mites and the total wombat population N on the centre-plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y = 0$$\end{document} , in the case in which this single dose of treatment has been applied, at the six different times indicated on the Figure and in its caption. Also shown in these two diagrams, for comparison, is the case when no treatment at all has been applied, so that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 0$$\end{document} throughout the entire environment. These untreated results are sketched with (purple) chain-dot lines.
Without any treatment, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r = 0$$\end{document} , Fig. 2 shows that the “plus” endemic state, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{eq} = 0.0324$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{eq} = 0.4663$$\end{document} , is the stable steady-state situation for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _F = 1/19$$\end{document} (the prevalence of observable mange in the wombat population for this case is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{H,eq} / N_{eq} = 1.2$$\end{document} percent). Thus, although the results were started from a perturbation to the mite-free equilibrium (16), the solution eventually makes a transition to the “plus” endemic state (24). The linearized solution of Sect. 4 is therefore not capable of following this transition between steady states (similar to the situation in Fig. 6), but the nonlinear method of Sect. 5 is able to cope with this strongly nonlinear effect, and the results are shown in Fig. 8 for both the completely untreated case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0$$\end{document} as well as the situation involving the single treatment over the time interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$200< t < 250$$\end{document} .
The use of a single treatment at the centre of the region evidently does not have a great effect on either mite numbers or overall wombat survivability. This is clear from Fig. 8, since after sufficient time has passed, there is almost no difference between the populations with treatment (blue, solid lines) and those without (purple, chain-dot lines). Nevertheless, at the second time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 240$$\end{document} shown in part (a), the treatment has clearly had a strong, short-term effect, since the peak in the mites for the untreated case is replaced by an inverted trough in the case with treatment. It is interesting to observe that, since the solution has been started from an initial condition near the unstable zero-mites equilibrium (16), mite numbers begin to grow exponentially at early times, just as predicted by the linearized solution in Sect. 4, before eventually undergoing nonlinear saturation at about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 0.2$$\end{document} at the time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 720$$\end{document} , and then making the transition to the “plus” endemic state for which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 0.0324$$\end{document} .
In an attempt to bring about a significant reduction in mite numbers F, repeated treatments over the same treatment region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_T$$\end{document} have been investigated, as detailed in (51). This has been achieved in this study by choosing the time-dependent component \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_T (t)$$\end{document} in (51) to be a switching function that activates at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=200$$\end{document} days and switches off after 50 days at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=250$$\end{document} , as in Fig. 8; however, it now reactivates every 250 days after this time. This creates the time-dependent switching function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_T (t)$$\end{document} illustrated in Fig. 9 over the total time period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0< t < 2400$$\end{document} considered here.Fig. 10. Comparison of treated case (solid line—blue online) with untreated case (chain-dot line—purple online). The parameters are the same as for Fig. 8. This figure illustrates: a the mite numbers F, and b the total wombat numbers N at times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=120$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 240$$\end{document} , 720, 1320, 1920 and 2400 days, on the centre-plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y=0$$\end{document} . The initial condition consisted of a perturbation to the mite-free state (16). Multiple treatment occurred over a square centred at the origin \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(x,y) = (0,0)$$\end{document} of the environment, lasting 50 days and repeating every 250 days. There was no ambient treatment, so that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A = 0$$\end{document} with special rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_T = 0.9/7$$\end{document} over the target region (color figure online)
The results of applying this repeated treatment are illustrated in Fig. 10. Here, the parameters are all the same as in Fig. 8, and the treatment region is chosen to be the same square patch \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< x, y < 1$$\end{document} located around the centre of the wombat environment. The ambient treatment rate is again \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A = 0$$\end{document} and the targeted treatment rate is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_T = 0.9/7$$\end{document} over the treatment patch. Now, however, that treatment is repeated periodically in time, as outlined in the schedule shown in Fig. 9. Accuracy was ensured by choosing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M,N = 121,121$$\end{document} Fourier coefficients in the numerical solution.
It is evident from Fig. 10 that the repeated treatment does indeed reduce the peak mite numbers over the approximate time interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$240< t < 1320$$\end{document} during which the solution is making its nonlinear transition from the mite-free state (16) to the “plus” endemic state (24). The repeated treatments have resulted in the formation of a “dimple” in the centre-line mite profile F(x, 0, t) in part (a), focussed at the origin. In addition, the overall wombat numbers in part (b) remain slightly above those for the untreated case. However, after sufficient time has passed, only slight advantages remain from the repeated treatment over the centre of the environment.Fig. 11. Some spatial patterns for selected times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=120$$\end{document} , 720 and 2400, for the same case illustrated in Fig. 10 with multiple treatment periods. This figure illustrates: a the mite numbers F, and b The total wombat numbers N for the treated case. There was no ambient treatment ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A = 0$$\end{document} ) and the treatment rate was \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_T = 0.9/7$$\end{document} over the central target region
Figure 11 illustrates further the way in which the repeated treatment over the rectangular target region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< x, y < 1$$\end{document} has resulted in a dent near the centre of the mite profile. Results are shown at the three times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 120$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=720$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=2400$$\end{document} days, and the indented mites profile is clearly evident at the second time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=720$$\end{document} in Fig. 11a. Evidently, the repeated application of treatment with the rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_T = 0.9/7 = 0.1286$$\end{document} (taken from Table 1) over the centrally-located square region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1< x, y < 1$$\end{document} has resulted in a substantial reduction in mite numbers over this special targeted portion of the environment for times about 720 days. Nevertheless, the long-term benefit of this treatment is fairly minor. This is a disappointing outcome, and suggests that an alternative treatment strategy should be investigated instead.
Variable treatment rate at a corner
Fig. 12. Some spatial patterns for selected times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=120$$\end{document} , 720 and 2400, for the same (treated) case illustrated in Fig. 10 with multiple treatment periods. There was no ambient treatment ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A = 0$$\end{document} ) and the treatment rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_T = 0.9/7$$\end{document} was applied over the square corner region indicated with (red) solid lines. This figure illustrates: a the mite numbers F, and b the total wombat numbers N (color figure online)
In this model, we have assumed that there is no movement across the boundaries of the spatial environment. For simplicity, we assumed a rectangular environment \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-L< x < L$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-B< y < B$$\end{document} , and the impermeability of these borders was expressed by the boundary conditions (7), (8). These boundaries might physically correspond to rivers, for example, that the wombats are not free to cross. This suggests that it might be possible to protect a group of wombats by moving the targeted treatment zone away from the centre, where the original group of mites was dropped at the initial time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 0$$\end{document} and where mite numbers are consequently greatest, and re-locating it to a corner of the environment.
Accordingly, we consider in Fig. 12 a situation in which the target region has been moved into one corner of the environment. In the notation of condition (52), we have chosen a target region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_T$$\end{document} with centre \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( x_0, y_0 \right) = (4,4)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_T = y_T = 1$$\end{document} ; this is a square \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3< x < 5$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3< y < 5$$\end{document} nestled into a corner of the environment, with two sides protected from wombat movement by the conditions (7), (8). To aid understanding, this region is sketched at the bottom of each diagram using a heavy solid line (red online). Figure 12a shows the mite numbers F over the entire environment at the three times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 120$$\end{document} , 720 and 2400 days, and the total wombat numbers N are illustrated in part (b) at those same times. These results were generated using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(M,N) = (121,121)$$\end{document} Fourier coefficients.
It can be seen, particularly from the last diagram at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=2400$$\end{document} in part (a), that mite numbers have been very substantially reduced over the target region, even although there is not much change over the rest of the environment. The total wombat numbers shown in part (b) also rise slightly in that region.Fig. 13. Some spatial patterns for selected times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=120$$\end{document} , 720 and 2400, for the same (treated) case illustrated in Fig. 10 with multiple treatment periods. There was no ambient treatment ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A = 0$$\end{document} ) and the treatment rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_T = 0.9/7$$\end{document} was applied over a larger square corner region, indicated with (red) solid lines. This figure illustrates: a the mite numbers F, and b the total wombat numbers N (color figure online)
Figure 13 continues this investigation of the effect of the special target region being located at a corner of the environment, with two sides prevented from the effects of wombat migration. Here, the situation is the same as for Fig. 12, except that the size of the corner target region has been increased to the square \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0< x < 5$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0< y < 5$$\end{document} which now occupies one quarter of the environment region. This square is sketched with solid lines (red online) at the bottom of each diagram. Repeated treatment is given over this region as in Fig. 12, following the schedule outlined in Fig. 9. As previously, there is no ambient treatment, so that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A = 0$$\end{document} , and the treatment rate is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_T = 0.9/7 = 0.1286$$\end{document} in the targeted region. To ensure high numerical accuracy, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(M,N) = (121,121)$$\end{document} Fourier coefficients are used in the spectral solution method (45). Sustained reductions in mite numbers, and corresponding increases in total wombat populations, are achievable in the protected target zone; by contrast, when the target zone was centrally located, it was not found possible to maintain the reductions in mite numbers that were achieved at early times, to the same extent.Fig. 14. Some spatial patterns for selected times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=120$$\end{document} , 720 and 2400, with multiple treatment periods applied in the same larger corner region as in Fig. 13. Here, however, there is an ambient treatment rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A = 0.01$$\end{document} across the entire region, with the treatment rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_T = 0.9/7$$\end{document} from Table 1 used in the target region (indicated with red solid lines). This figure illustrates: a the mite numbers F, and b the total wombat numbers N (color figure online)
Finally, we consider the same case as in Fig. 13, except that, now, some ambient treatment over the entire region is supposed to be present, at the low rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_A = 0.01$$\end{document} . The results are displayed in Fig. 14. From Figs. 2 and 3 , the “plus” endemic state is still the stable steady state, and this has been confirmed from the linearized solution in Sect. 4. This equilibrium has fomite number density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F^{(+)}_{eq} = 0.0060$$\end{document} and total wombat number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^{(+)}_{eq} = 0.7622$$\end{document} , as calculated from (18) and (19). As previously, the calculations are started from the zero-mites steady state (16) for which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{eq} = 0.8$$\end{document} , perturbed by the injection of a patch of mites of population size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon = 0.2$$\end{document} . The solutions in Fig. 14 were computed with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(M,N) = (121,121)$$\end{document} Fourier coefficients.
From the last image at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 2400$$\end{document} in Fig. 14a, it is clear that the primary aim of protecting the wombats in the corner target region (indicated on each diagram) has essentially been accomplished. The mite number density F remains very low over that region, over a sustained period of time. Furthermore, from Fig. 14b, the total wombat population N continues to remain closer to the mite-free equilibrium value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{eq} = 0.8$$\end{document} , rather than falling to the endemic equilibrium value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^{(+)}_{eq} = 0.7622$$\end{document} . This is a potentially encouraging outcome for wombat survival.
Discussion
Sarcoptes scabiei causes a significant disease burden, on both humanity (World Health Organisation 2020) and a variety of mammalian species across the globe (Astorga et al. 2018). Research addressing management of mange is therefore of value for mitigating this burden. Here, we have sought to explore a bare-nosed wombat population living with endemic mange, and how that population responds to a standard treatment course (Martin et al. 2019) in a restricted treatment area.
The model system we have investigated here predicts three biologically meaningful steady-state configurations. The first is a total extinction state (15), in which the populations of both the mites and the wombats all become zero. This would represent a disastrous outcome for wombat conservation in that area, and must be avoided. Fortunately, we have shown in Sect. 4 that this extinction steady-state possibility is unstable, so that, even if wombat populations become gravely low, the dynamical behaviour of this system will not inevitably drive them to extinction. The second possible steady state is a mite-free configuration. This is surely the most desirable outcome, since the mites are all eliminated and only the wombats remain. An important consideration, then, is to determine under what conditions this mite-free state is stable, which from a practical point of view would imply that if mites were to be introduced, they would eventually all die out. In general this important question can only be answered numerically; however, when the two death rates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _L$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _H$$\end{document} are equal, we have managed to derive, in Theorem 1, a condition which is both necessary and sufficient for the stability of this most desirable steady state. That condition (39) is illustrated graphically in Fig. 3. There is also a third endemic steady state in which both wombats and mites coexist. This is clearly a less desirable outcome, since wombats would suffer the ill effects of mange. (Mathematically, there is also a fourth steady state, but it involves negative populations, and so makes no biological sense).
A linearized approximation to these governing equations was derived in Sect. 4, and the behaviour of its solutions was analyzed. It provides a check on the reliability of computer codes that are devised to solve the governing equations, and it is also used to assess the stability of the equilibrium states. Importantly, there are situations in which small changes to background conditions, as expressed through the parameters in the model, can bring about large-scale qualitative changes to the outcomes for the wombats. In particular, Figs. 2 and 3 indicate that mite-free steady states are possible with appropriate treatment rates, if the treatment can be made to cover the entire environment. By contrast, results such as those in Fig. 11 make clear the fact that, in general, treatment is unlikely to be greatly effective when it only reaches some fraction of the wombat habitat. In practice, it is difficult in the field to know how successful disease management has been, since it is usually difficult to locate all the wombat burrows, or to be certain how efficiently the burrow flaps are working, and so on. At present, there are limited robust field studies that document the management of S. scabiei in wombat populations, outside of the study undertaken by Martin et al. (2019).
The linearized solution of Sect. 4 confirmed the veracity of the non-linear results, in circumstances where a population remained reasonably close to an assumed equilibrium state. However, the linearized solution is not capable of predicting a transition from one equilibrium state to another, since it is predicated on the idea that populations consist of a small perturbation to a given equilibrium situation. Under such circumstances, linearization fails after sufficient time has elapsed, since linearized solutions can either only decay or grow exponentially. Nevertheless, they are reasonably accurate even in unstable circumstances at early times, and our linearized solution in Sect. 4 was able to confirm our nonlinear result in unstable situations such as that illustrated in Fig. 6, too, but only at very early times. Later, however, the linearized solution continues to increase exponentially in size, so violating the conditions under which linearization itself is valid. The nonlinear solution, on the other hand, retains its validity at these later times; the solution becomes saturated by nonlinear effects, and at least in some instances, it transitions toward a different equilibrium state.
The behaviour of the nonlinear solution with time can be very difficult to predict, even if it ultimately tends to one of the known equilibrium states in Sect. 3. Its time-dependent behaviour is often sensitively dependent on initial conditions, and this has been illustrated here. The nonlinear solution can undergo large-amplitude excursions, with explosive growth of mite numbers in the shorter term. This was seen in the test case illustrated in Fig. 6, and confirmed with parameters relevant to wombat–mite interactions in Figs. 10, 11, 12 and 13. This could have important ramifications for how biological systems react to sudden changes in their environment, for example. In some circumstances it may be a vital consideration, and is deserving of careful attention in systems such as that studied here; nevertheless, the inherent nonlinearity of the system makes such behaviour very difficult to quantify.
In cases where the naturally-occurring (stable) situation is one in which both wombats and mites co-exist endemically in the environment, even with a low level of ambient treatment against mites, we have investigated situations in which additional higher levels of treatment are applied, but only to certain limited regions in space. The aim, of course, is to reduce mite numbers and boost wombat survivability, at least in those special targeted regions. We have demonstrated, however, that while there is indeed some shorter-term reduction in mite number density in those target regions, the overall spatio-temporal behaviour of the system is somewhat more subtle and nuanced than perhaps expected.
When the targeted treatment region is located centrally within the mite-infested region, the benefits of aggressive treatment within that region, even over repeated time intervals, are rather mixed. Although there is an initial reduction in mite numbers, this becomes considerably less pronounced at later times; furthermore, the increase in overall wombat numbers is almost negligible. However, our calculations have demonstrated that, when it is possible to take advantage of naturally-occurring borders or boundaries within the environment, such as rivers that the wombats are not able to cross, then considerably better outcomes may be possible by locating the targeted treatment region judiciously next to these features. In that case, we have shown that it is possible to achieve almost zero mite numbers within those target regions, and particularly when even a quite low level of ambient treatment is possible throughout the whole region. Furthermore, it is possible to achieve long-term wombat numbers that remain slightly above what would be expected from the equilibrium analysis alone (as in Sect. 3), and this was illustrated in Fig. 14. Our model has assumed the wombat habitat to be a perfect rectangle, and although this is highly unlikely to be true in the field, the local reductions in mite numbers we were able to achieve in Fig. 14, near one corner of the region, should be able to be observed in practice. This is because our result is local to a corner and does not rely on changes far away, so that a naturally-occurring similar corner with real impermeable boundaries offers the likelihood of major benefits to wombats in that region.
In this paper, we have modelled wombat (and burrow) distribution as relatively uniform across our simulated spatial environment, so as to limit variables in our model and to support tractable mathematical analysis of the problem. However, wombats are known to have habitat preferences in where they distribute their burrows, and this could be a topic for further research, perhaps even using individual based models. Our field research suggests that wombats generally have fairly fixed and stable home ranges, so our model is likely to represent a reasonable approximation.
There is clearly scope for future work to investigate further strategies to protect wombat populations, based on the location and timing of the treatment regime. Ideally, the entire problem might even be thought of as an optimal control system that is aimed at maximizing wombat numbers, using the treatment function r(x, y, t) as the control variable. This would, however, be an enormously ambitious theoretical undertaking that would also make high demands on computer resources, rendering it currently infeasible. Nevertheless, with the careful use of a simulation tool similar to that developed in this paper, combined with detailed knowledge of the terrain and the location of wombats within it provided by dedicated operators in the field, long-term protection of wombat populations from the most severe effects of infestation by sarcoptic mange may be achievable.
Supplementary information
Highly efficient vectorized code has been developed for the implementation of the nonlinear algorithm presented in Sect. 5 and is presented in Supplementary Material. Further numerical details concerning this code are given in the “Appendix A”.
Supplementary Information
Below is the link to the electronic supplementary material.Supplementary file 1 (pdf 98 KB)
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Arlian LG Morgan MSA review of Sarcoptes scabiei: past, present and future Parasites Vectors 201710.1186/s 13071-017-2234-1PMC 547775928633664 · doi ↗ · pubmed ↗
- 2Arlian LG Runyan RA Achar S Estes SA Survival and infestivity of Sarcoptes scabiei var. canis and var. hominis J Am Acad Dermatol 19841121021510.1016/S 0190-9622(84)70151-46434601 · doi ↗ · pubmed ↗
- 3Astorga F, Carver S, Almberg ES, Sousa GR, Wingfield K, Niedringhaus KD, Van Wick P, Rossi L, Xie Y, Cross P, Angelone S, Gortázar C, Escobar LE (2018) International meeting on sarcoptic mange in wildlife, June 2018, Blacksburg, Virginia, USA. Parasites and Vectors, vol 11, p 449. 10.1186/s 13071-018-3015-110.1186/s 13071-018-3015-1PMC 609081330075742 · doi ↗ · pubmed ↗
- 4Banks SC Skerratt LF Taylor AC Female dispersal and relatedness structure in common wombats (Vombatus ursinus)J Zool Lond 200225638939910.1017/s 0952836902000432 · doi ↗
- 5Beeton NJ Carver S Forbes LKA model for the treatment of environmentally transmitted sarcoptic mange in bare-nosed wombats (Vombatus ursinus)J Theor Biol 201946246647410.1016/j.jtbi.2018.11.03330502410 · doi ↗ · pubmed ↗
- 6Browne E Driessen MM Ross R Roach M Carver S Environmental suitability of bare-nosed wombat burrows for Sarcoptes scabiei Int J Parasitol Parasites Wildl 202116374710.1016/j.ijppaw.2021.08.00334434693 PMC 8374697 · doi ↗ · pubmed ↗
- 7Death CE Taggart DA Williams DB Milne R Schultz DJ Holyoake C Warren KS Pharmacokinetics of moxidectin in the southern hairy-nosed wombat (Lasiorhinus latifrons)J Wildl Dis 20114764364910.7589/0090-3558-47.3.64321719829 · doi ↗ · pubmed ↗
- 8Driessen MM Dewar E Carver S Gales R Conservation status of common wombats in Tasmania I: incidence of mange and its significance Pac Conserv Biol 20212810311410.1071/PC 21007 · doi ↗