A system dynamics model of the population dynamics of Oestrus sp. (Diptera: Oestridae) infesting Iberian ibex, Capra pyrenaica

Abstract

We modelled the population biology of Oestrus sp. parasitising the Iberian ibex (Capra pyrenaica), using a system dynamics approach. Levels included in the model were as follows: first-instar larvae, third-instar larvae, pupae and resilient pupae. The values used for flows and auxiliary variables were obtained from the scientific literature, including research on both Oestrus ovis and O. sp.; we assumed that these two species are closely related, both morphologically and from a biological point of view. Simulations provided by our model allowed us to estimate the minimum monthly production of first-instar larvae and thus to establish periodic population dynamics, together with predictions for the number of larval generations. The model predicted 3–4 larval generations per year. The model proved to be very sensitive to minimal changes in a number of variables, especially the first-instar larval production rate. Despite its limitations, this methodology could be a versatile tool for studying the population dynamics of this kind of parasites, and for simulating the effects of control programmes.

Keywords:

• Population biology
• modelling
• system dynamics
• Oestrus sp.
• Oestridae

Introduction

Oestrids (Diptera: Oestridae) are obligate parasites that remain for weeks or months in the naso-pharyngeal tract, frontal and horn sinuses (Oestrus spp.) (Figure 1), stomach (Gasterophilus spp.), inner organs or subcutaneous tissues (Hypoderma spp.) of their respective hosts. They cause important economic loses in livestock in both developing and developed countries (Otranto & Stevens 2002), and under certain conditions can also infect humans (Anderson 2006b).

Figure 1. The life cycle of Oestrus sp.

Beyond the mechanical damage caused by Oestrus ovis Linnaeus, 1758 larvae in host mucosae, other alterations may appear, including dilatations, degenerative changes and squamous metaplasia in nasal passages, as well as leukocyte infiltration – mainly eosinophils – or submucosal edema (Jagannath et al. 1989). Chronic sinusitis, probably due to secondary infections associated to these myasis, has often been reported in bighorn sheep (Ovis canadensis Shaw, 1804) and may affect host crania and even cause death (Allen & Bunch 1982).

Oestrus ovis parasitises livestock (mainly sheep) the world over and has also been reported in bighorn sheep, European mouflon (Ovis aries Linnaeus, 1758), argali (Ovis ammon (Linnaeus, 1758)), Alpine ibex (Capra ibex ibex Linnaeus, 1758), Asiatic ibex (Capra sibirica (Pallas, 1776)) and white-tailed deer (Odocoileus virginianus (Zimmermann, 1780)) (Grunin 1957; Capelle 1966; Wetzel & Bauristhene 1970; Moreno et al. 1999; Colwell et al. 2006). In Eurasia Oestrus caucasicus Grunin, 1948 is known to parasitise the Kuban tur (Capra caucasica Güldenstaedt and Pallas, 1783), the Asiatic ibex, and also was reported infesting the Iberian ibex (Capra pyrenaica Schinz, 1838) (Grunin 1957; Minar et al. 1985; Pérez et al. 1996). Recent molecular studies based on Cytochrome Oxidase I (COI) and 28S rDNA sequences confirmed that the species infesting Iberian ibex is other than O. ovis (Moreno et al. in press). Nevertheless, since we have not yet obtained and genetically characterised larvae from Kuban tur, Daghestan tur (Capra cylindricornis (Blyth, 1841)) and Asiatic ibex, we cannot conclude that the species infesting Iberian ibex is O. caucasicus or, even, a new species. For this reason, hereafter, we call it Oestrus sp. Regarding morphological and biological features, this latter parasite species is very similar to O. ovis (Grunin 1957; Zumpt 1965; Colwell 2001; Guitton et al. 2001; Pérez et al. 2006).

Monitoring populations of these parasites when working on infected wild hosts is a difficult, time-consuming task. Several non-invasive methods are available for diagnosing oestrosis in wild hosts, such as Enzyme-Linked Immunosorbent Assay (ELISA) (Arias et al. 2014), pharyngeal endoscopy (Oksanen et al. 1992) or computed tomography (Fidalgo et al. 2015). Nevertheless, for epidemiological studies and a proper identification of the parasite species involved, larvae are usually obtained after necropsy (invasive method). In addition, pupae develop in the substrate at a certain depth and the behaviour of adult flies makes it extremely difficult to observe them in vivo, taking into account that the life cycle of the species involved in our study occurs in high mountain habitats. Therefore, our knowledge of the life cycle and demographic patterns of oestrids infesting wild hosts is scarce and fragmented, to the point that most of their demographic parameters are still unknown.

System dynamics (SD) modelling is a methodology used for simulating and analysing the behaviour of complex models (Aracil & Gordillo 1997). SD acts as an intermediate stage between reality and its mathematical representation and can be used to predict by means of simulations important temporal variations (Martín-García 2003). From a mathematical point of view, a SD model consists of a set of first-order ordinary differential equations. Under the circumstances described above, models or simulations may help to obtain a structured vision of the infestation process and its critical aspects, and provide possible approaches for managing these infections. Within this context, the basic goal of a model is to understand how the different variables involved keep the system in a balanced stationary state, and how the system will react if a variable is modified.

There are some antecedents in the literature regarding the modelling population dynamics of Diptera using data obtained experimentally, although most refer to the Calliphoridae and Muscidae families. Some of these studies have been focused on constructing life tables that allow important demographic parameters such as the net reproductive rate (R₀), the generation time and the instantaneous growth rate to be derived, and to predict the future structure of the studied population (Chi & Liu 1985; Chi 1988, 1997; Gabre et al. 2005).

Given that insect developmental and survival rates are determined by environmental factors such as temperature, it is not surprising that most insect population models try to quantitatively extrapolate a rate of population increase using available climate data (Atzeni et al. 1994; Wall et al. 2000). However, it must be taken into account that these models usually include the calculation of degree-days under the assumption that above a minimum threshold (basal temperature) developmental rates are directly proportional to temperature, and that beyond this threshold development does not occur (Moon 1983; Preuss KP 1983).

Nevertheless, the control of the common green bottle fly, Lucilia sericata (Meigen, 1826), has been approached with the help of both deterministic (assuming no inter-individual variability) and stochastic models, which assume a certain degree of variability between individuals and cohorts (Wall et al. 1992, 1993a,b; French & Morgan 1996; Fenton & Wall 1997; Fenton et al. 1997, 1998a,b). Within the context of integrated insect-pest management, the monitoring of insect populations has also been performed by using mixed generalised lineal models (Candy 2000).

Little life cycle modelling has been carried out in oestrid research. In order to predict critical periods with the greatest numbers of attacks on hosts (larvipositions) and to design preventive and control strategies, the degree-days and the monitoring of the prevalence and intensity of oestrosis in sheep caused by O. ovis from Namibia have been used to address the question of development and/or survival of different parasite stages (adult, larva and pupa) (Biggs et al. 1998). On the other hand, recent work has reported a life table for an O. ovis population and predicted a rate of population growth of 1.25 per generation, under natural conditions (Cepeda-Palacios et al. 2011).

In our study, the life cycle of Oestrus sp. infesting Iberian ibex (Capra pyrenaica) is modelled using a SD approach. Given the similarity of this species to O. ovis, we integrated available information for both species from larvae (parasitic phase), pupae and adult flies (free-living phases). The goal of our study is to estimate the minimum number of first-instar larvae produced by gravid females which is needed to maintain stable (periodic) population dynamics, assuming high rates of larval and pupal mortality. This information, after being validated empirically, could become crucial within control programmes for this parasitosis.

Material and methods

Elements and structure of the system dynamics model

In SD, a flow diagram (or Forrester’s diagram) consists of three main elements: (i) levels, the state variables of the system with their respective values in each situation as the difference between input and output flows; (ii) flows, used to establish the relationship between the different levels; and (iii) auxiliary variables, which integrate information channels between levels and flows and help to define the flow variables, thereby making the model more manageable and easier to understand. Our model includes four levels: L1 (first-instar larvae), L3 (third-instar larvae), pupae and resilient pupae; eight flows: new L1, L1 death, from L1 to L3, L3 death, from L3 to pupae, pupae death, resilience and unmeeting; and 20 auxiliary variables: larviposition, L1 death rate, L1 number, delay 1, control cycle L1 to L3, delay control cycle L1 to L3, time L1 to L3, L3 death rate, L3 number, control cycle L3 to pupae, delay 2, delay control cycle L3 to pupae, time L3 to pupae, pupae number, pupae death rate, delay 3, resilient time, control cycle pupae to resilient pupae, unmeeting rate, and offspring potential. From these auxiliary variables, larviposition (that is, the rate at which gravid females laid first-instar larvae within the host nostrils) is one of the most important with regards to the population dynamics of this fly, taking into account the very short life expectancy of adults.

Values assigned to the relationships between levels and auxiliary variables of the model were obtained from the literature. Adult life expectancy was considered as the mean obtained for males and females: 0.31 ± 0.13 months (Pérez et al. 2006). The biotic potential for Oestrus sp. gravid females was considered to be the same as previously described for O. ovis (500 larvae), but, taking into account that only 30–50 larvae are laid in each larviposition (Cobbett & Mitchell 1941; Zumpt 1965), we included the auxiliary variable “Larviposition”, which was parameterised after making the model periodic. In addition, we have taken into account two larval instars (L1 and L3) in order to simplify the system. Assuming overall larval mortality of 90% (Rogers & Knapp 1973) and depending on climatic conditions (mainly temperature), the complete larval phase (L1–L3) lasts from 20 days to a number of months (Dorchies et al. 1999; Tabouret et al. 2001). Thus, even if it is included, the intermediate L2 instar does not become informative or significant to the model. In the model, we assumed a long larval period of 25 days (0.83 months): “Time L1 to L3”, which is also an auxiliary variable, as the arithmetic mean of the values obtained by Dorchies et al. (1999) and Tabouret et al. (2001). In addition, in vitro culture experiments showed that Oestrus sp. post-feeding third-instar larvae have a mortality rate of 10.6% and that they take 1.28 days (“Time L3 to Pupae”) on average to pupate, the pupal stage then lasting for 30 days (Pérez et al. 2006); all of these values were included in the model.

The emergence rate of adult flies reaches 43.6% (assuming a mortality rate of 56.4% during the pupal phase) with a sex ratio not differing significantly from 1:1 (Pérez et al. 2006). However, we assumed that 17.9% of adult flies emerged with vestigial wings and would probably not reproduce (Pérez et al. 2006). Because we have no information about the natural mortality of adults (e.g. by predation) nor regarding the success of mating and the rate of larviposition by gravid females (e.g. the success regarding different larviposition attempts), we have integrated the mortality rate of the pupal phase, the sex ratio, and the rate of vestigial adults into the level “resilient pupae”. The offspring potential reaches a value of 250 (since only about the 50% of adult flies are females).

In order to emulate the seasonal fluctuations in the life cycle, particularly those regarding the number of larvae belonging to each instar (Pérez et al. 1996), we have introduced different fixed delays in the model, just prior to passing to the next level, and three control cycles as well. These control cycles are defined as step functions and are also important auxiliary variables of the model.

The time unit used in this work was the month. The model of the life cycle for Oestrus sp. was developed with the software Vensim PLE 6.2 (Ventana Systems, Inc. 2003). In the first simulation, using an initial value of eight L1, 18 L3, 16 pupae and eight resilient pupae, we looked for the conditions that guaranteed the periodic behaviour of the population – that is, a rate of increase over time equal to or very near to zero (Caughley & Gunn 1996). These initial values are based on data on intensity of parasitation described by Pérez et al. (1996). The time length for the simulation was 25 years, for which we used the Runge Kutta method for numerical integration, with a step of h = 0.015625.

Finaly, after estimating the value of larviposition which makes the model stable, we carried out a sensitivity analysis of the model based on 200 simulations using a random–uniform probability distribution for larviposition rate values between 0.25 and 12.5 larvae/gravid female/month.

Results

The general structure of the model implemented with Vensim PLE 6.2 is represented by the equations of the model (Table I) and depicted in Figure 2.

Table I. Equations describing the model. L1: first-instar larvae; L3: third-instar larvae; INTEG: integral.

[1]	Pupae number = pupae/75
[2]	L1 number = L1/82
[3]	L3 number = L3/75
[4]	Offspring potential = resilient Pupae × 250 × control cycle L1 to L3
[5]	Resilient pupae = INTEG (resilience-unmeeting, 8)
[6]	Control cycle pupae to resilient pupae = PULSE TRAIN(9, 1,12,500)
[7]	From L3 to pupae = delay control cycle L3 to pupae × delay 2/time L3 to pupae
[8]	Resilience = control cycle pupae to resilient pupae × delay 3/resilient time
[9]	Delay control cycle L3 to pupae = DELAY3(control cycle L3 to pupae, 0.1)
[10]	Control cycle L3 to pupae = 0.03 × PULSE TRAIN(0,1,12,500)+0.005 × PULSE TRAIN(1,2,12,500)+0.0075 × PULSE TRAIN(2,1,12,500)+ 0.02 × PULSE TRAIN(3,1,12,500)+0.05 × PULSE TRAIN(4,1,12,500)+0.03 × PULSE TRAIN(5,1,12,500) +0.015 × PULSE TRAIN(6,1,12,500)+0.055 × PULSE TRAIN(7,1,12,500)- 0.04 × PULSE TRAIN(8,1,12,500)+0.075 × PULSE TRAIN(9, 1,12,500)−0.05 × PULSE TRAIN(10, 1,12,500)
[11]	new L1 = offspring potential × larviposition
[12]	Control cycle L1 to L3 = 0.25 × PULSE TRAIN(1,1,12,500)−0.25 × PULSE TRAIN(2,1,12,500)+ 0.3 × PULSE TRAIN(3,1,12,500)−0.9 × PULSE TRAIN(4,1,12,500)+1.2 × PULSE TRAIN(5,1,12,500)−0.2 × PULSE TRAIN(6,1,12,500)−0.5 × PULSE TRAIN(7,1,12,500)+0.6 × PULSE TRAIN(8,2,12,500)+0.8 × PULSE TRAIN(11, 1,12,500)
[13]	L1 death = L1 × L1 death rate
[14]	Resilient time = 1
[15]	From L1 to L3 = delay 1 × delay control cycle L1 to L3/time L1 to L3
[16]	Time L1 to L3 = 0.83
[17]	Time L3 to pupae = 0.043
[18]	Delay control cycle L1 to L3 = DELAY3(control cycle L1 to L3, 1)
[19]	L1 = INTEG (+new L1-L1 death-from L1 to L3 + 60.8)
[20]	L3 = INTEG (from L1 to L3-L3 death-from L3 to pupae+400.18)
[21]	Larviposition = 0.02
[22]	Delay 2 = DELAY3(L3, 0.5)
[23]	Delay 1 = DELAY3(L1, 0.8)
[24]	Delay 3 = DELAY3(pupae, 0.25)
[25]	Unmeeting rate = 0.4
[26]	Unmeeting = resilient pupae × unmeeting rate
[27]	L1 death rate = 0.9
[28]	Pupae = INTEG (from L3 to pupae-pupae death-resilience+5.16)
[29]	Pupae death = pupae × pupae death rate
[30]	Pupae death rate = 0.564
[31]	L3 death = L3 × L3 death rate
[32]	L3 death rate = 0.106

Figure 2. A Forrester’s diagram representing the elements of the model together with their relationships. L1: first-instar larvae; L3: third-instar larvae.

Each of the level variables represents a sector of the population and has an input flow (births) and two output flows (deaths); the velocities or rates at which these flows operate are density dependent. These rates are regulated by the control variable “delay control cycle L1 to L3” to pupae” for the period L3 pupa.

Starting from a given initial population as described below, we obtain a stable or stationary population (system) dynamics when the mean monthly rate reaches a value of five larvae per adult gravid female and month. Moreover, the model proved to be very sensitive to small changes in the value of different auxiliary variables, especially that for larviposition (Figure 3).

Figure 3. The sensitivity of the model, obtained after 200 simulations (larviposition ranging from 0.25 to 12.5 larvae/gravid female/month) regarding L1 (first-instar larvae) numbers (y-axis). On the x-axis, 24 and 36 represent December, and 25 to 35 represent January to November, respectively.

The model predicts 3–4 larval generations per year (Figures 4–5). In general, the values of both L1 and L3 predicted by the model are very similar to those described by Pérez et al. (1996). The model has a cyclic behaviour and, then, the predicted dynamics of the population is periodic (Figure 6).

Figure 4. Dynamics of L1 (first-instar larvae) numbers. In red (bottom legend), number of L1 observed (Pérez et al. 1996). In blue (upper legend), number of L1 predicted by the model. In the x-axis, 24 and 36 represent December, and 25 to 35 represent January to November, respectively.

Figure 5. Dynamics of L3 (third-instar larvae) numbers. In red (bottom legend), number of L3 observed (Pérez et al. 1996). In blue (upper legend), number of L3 predicted by the model. In the x-axis, 24 and 36 represent December, and 25 to 35 represent January to November, respectively.

Figure 6. Limit cycle of the model. As the curve is closed, the long-term behaviour of the model is periodic. L1: first-instar larvae; L3: third-instar larvae.

Discussion

Monitoring oestrosis throughout the entire life cycle of these parasites would be extremely difficult. Thus, studies usually focus on larval-parasitic (epidemiology) or adult phases in the laboratory (after larval culture), or study adult biology in the field (Anderson 2006a). However, simulation gives us an opportunity to relate all the different phases of an oestrid’s life cycle in order to obtain predictive values under specific conditions or values for each parameter.

Matrix population models (Caswell 2001) have many applications in population dynamics. They basically operate with a matrix including the age population structure for each sex class and two parameters: one including the survival rate and the other the reproductive rate for each sex and age class. Nevertheless, in our case: (i) we don’t know the sex of the larvae; (ii) they do not participate in reproduction (therefore, the matrix would contain many zeros); and (iii) we don’t know the mating success, the reproductive rate or the survival rate of adult flies.

SD allows us to better understand the system under study and the way in which its different elements are related each other. Validation of models, when empirical data are available, can be achieved easily. SD permits obtaining the value of one of its variables or elements indirectly. Moreover, there is specific software available for working with SD.

The seasonal pattern of prevalence and intensity of parasitation by Oestrus sp. larvae in Iberian ibex from the southern Iberian Peninsula included three generations of first-instar larvae per year (May, July and October), followed by successive generations of L2 and L3 (Pérez et al. 1996). In this paper, the highest values of mean intensity of parasitisation were achieved during the months December–January, and this could explain the “presumable” fourth peak in L1 production observed in these months, since it is not probable that adults would be active in a high mountain habitat during this period. The prediction made by our model (Figure 4) agrees with results obtained by Pérez et al. (1996).

On the other hand, under SD we cannot simulate levels that reach a value of zero, since this will cause the system to collapse. In real life in a high mountain habitat, we would not expect to find adult flies in the period between autumn and the following spring. During this period, the parasite population is maintained by larvae within their respective hosts and by pupae at a certain profundity in the substrate, both phases showing a certain degree of diapause. In the Ebro Valley (northeast Spain), Oestrus ovis adults were present and active between May and November, and a larval diapause from October to February was observed (Gracia et al. 2006). This is a subject that must be studied in greater detail.

Other important factors that are difficult to model because they are not fully understood are related to parasite embryonic development, its reproductive success (generally below its biotic potential; Colwell 2001), the behaviour of adult flies during reproduction, host seeking and larviposition related to the defensive behaviour of hosts (Anderson 2006a), the immune response of hosts (Arias et al. 2014), the influence of the programmes for treating and controlling oestrosis (Colwell et al. 2009), and even the number of hosts present in a given area. All of these factors will determine in each case the parasite carrying capacity and the density or abundance of adult flies in a given area. In other words, basic epidemiological variables such as prevalence and incidence (including reinfestations) also need to be included in models.

Therefore, we consider that this is only a first step towards predicting certain demographic values for parasites, which would otherwise be very difficult to obtain in the field or laboratory. Nevertheless, the results obtained in this work – e.g. the minimum larval production rate (five larvae per adult female and per month) required to maintain population periodicity – suggest interesting hypotheses that need to be tested experimentally. This kind of information could become essential for oestrosis control programmes.

Future research could be focused on each phase of the parasite’s cycle in order to discover how factors affecting both host and parasite – and internal and external to the system as a whole – operate and could generate more precise predictions. On the other hand, we should also consider the potential of stochastic models, which are more realistic but generate more uncertainty when predictions are made in the mid- or long term.

Aside from the host species involved, the geographical location and the respective climatological variables play an important role in the epidemiology of oestrosis. Other exogenous factors not taken into account in the model include the programmes developed for the treatment and control of these myiases in different parts of the world (Louw 1989; Pandey 1989).

Acknowledgements

The research activities of the authors are partially supported by the Junta de Andalucía: RNM118 research group. Comments by Dr. D. Colwell allowed us to greatly improve a preliminary version of our manuscript. The authors are also indebted to Michael Lockwood for his revision of the English style.