Nina Alphey
email: nina.alphey@zoo.ox.ac.uk

Michael B Bonsall
email: michael.bonsall@zoo.ox.ac.uk

Luke Alphey
email: luke.alphey@oxitec.com

Modeling resistance to genetic control of insects

Keywords:

The sterile insect technique is an area-wide pest control method that reduces pest populations by releasing mass-reared sterile insects which compete for mates with wild insects. Modern molecular tools have created possibilities for improving and extending the Modeling resistance to genetic control of insects page 3

sterile insect technique. As with any new insect control method, questions arise about potential resistance. Genetic RIDL ® technology (Release of Insects carrying a Dominant Lethal) is a proposed modification of the technique, releasing insects that are homozygous for a repressible dominant lethal genetic construct rather than being sterilized by irradiation. Hypothetical resistance to the lethal mechanism is a potential threat to RIDL strategies' effectiveness. Using population genetic and population dynamic models, we assess the circumstances under which monogenic biochemically-based resistance could have a significant impact on the effectiveness of releases for population control. We assume that released insects would be homozygous susceptible to the lethal genetic construct and therefore releases would have a built-in element of resistance dilution. We find that this effect could prevent or limit the spread of resistance to RIDL constructs; the outcomes are subject to competing selective forces deriving from the fitness properties of resistance and the release ratio. Resistance that is spreading and capable of having a significant detrimental impact on population reduction is identifiable, signaling in advance a need for mitigating action.

FOOTNOTE

RIDL ® is a registered trademark of Oxitec Limited, UK.

Introduction

Pest insects cause enormous harm to public health and economic damage to agriculture.

Arthropods act as vectors for many human and veterinary diseases and annually destroy around 18% of the world crop production [START_REF] Nicholson | Fighting the global pest problem: Preface to the special Toxicon issue on insecticidal toxins and their potential for insect pest control[END_REF]. The Sterile Insect Technique (SIT) is a pest population control method in which insects are mass-reared, sterilized by irradiation and released [START_REF] Dyck | Sterile Insect Technique: Principles and practice in Area-Wide Integrated Pest Management[END_REF][START_REF] Knipling | Possibilities of insect control or eradication through use of sexually sterile males[END_REF]. Some wild insects mate with the released sterile insects and produce fewer progeny as a result. If sufficient insects are released for a sufficient period, populations of pest insects can be suppressed or perhaps eliminated locally. The SIT has been successful against a range of target species, mainly agricultural pests [START_REF] Dyck | Sterile Insect Technique: Principles and practice in Area-Wide Integrated Pest Management[END_REF].

Genetic technology, such as RIDL ® (Release of Insects carrying a Dominant Lethal), could modify the SIT by releasing insects that are not irradiated, but instead are homozygous for one or more dominant lethal genetic constructs that are repressible during mass-rearing (Alphey et al., 2007a;[START_REF] Alphey | Engineering insects for the Sterile Insect Technique[END_REF][START_REF] Thomas | Insect population control using a dominant, repressible, lethal genetic system[END_REF]. This approach is also known as autocidal biological control [START_REF] Fryxell | Autocidal biological control: a general strategy for insect control based on genetic transformation with a highly conserved gene[END_REF]. Among the potential advantages claimed for genetic control over radiation-based SIT is the ability to engineer the timing or nature of the lethality. The preferred time of death varies between pests and contexts. Very early acting lethality would be preferable for pests where the larval stage causes the most damage; an embryonic lethal system has been developed in a model organism and work in a pest fruit fly is ongoing [START_REF] Horn | A transgene-based, embryo-specific lethality system for insect pest management[END_REF][START_REF] Schetelig | Development of an embryonic lethality system in Mediterranean fruit fly Ceratitis capitata[END_REF]. Late-acting lethality that delays the time of death until after a density-dependent mortality phase, such as has been developed for RIDL technology in mosquitoes, can greatly enhance the SIT's effectiveness against a population that is limited by density-dependent processes [START_REF] Atkinson | Analyzing the control of mosquito-borne diseases by a dominant lethal genetic system[END_REF][START_REF] Phuc | Late-acting dominant lethal genetic systems and mosquito control[END_REF]. Female-specific lethality would be advantageous for several reasons connected with facilitating separation of the sexes [START_REF] Dyck | Sterile Insect Technique: Principles and practice in Area-Wide Integrated Pest Management[END_REF]. Release of only males, rather than both males and females, has been shown to make the SIT much more effective [START_REF] Rendón | Medfly (Diptera:Tephritidae) genetic sexing: large-scale field comparison of males-only and bisexual sterile fly releases in Guatemala[END_REF]. The efficiency of releases is improved because the released males have to seek out wild females; if both sexes were released there would be potential for assortative matings among the released sterile males and females and this would reduce the effect on control of the natural population. For insect vectors where females transmit disease, such as mosquitoes, release of females would be highly undesirable. However, in the absence of genetic sexing mechanisms, the biology of most species does not permit sex separation on the large scale required for the SIT and male-only releases are currently the exception rather than the norm.

A female-lethal version of RIDL, with insects homozygous for one or more female-specific dominant lethal genetic constructs, has been constructed in several species [START_REF] Fu | Female-specific insect lethality engineered using alternative splicing[END_REF]Stainton et al., unpublished data). F 1 progeny of RIDL males and wild females inherit a dominant female-specific lethal gene; the F 1 females die, thereby reducing the reproductive potential of the wild population, but the F 1 males are viable and fertile. This provides a genetic sexing mechanism facilitating male only release, either by employing the femalelethal version of RIDL and withdrawing the repressor from the generation prior to release, or by combining a bisex-lethal system with female-lethality (with an independent means of repressing or inducing lethality) to permit male only release of bisex-lethal strains designed to kill progeny of both sexes in the field.

Insects are well-known for developing resistance to insecticides. Insecticide resistance has been reported in most major insect disease vectors and against every class of chemical insecticide and insecticidal crops [START_REF] Nicholson | Fighting the global pest problem: Preface to the special Toxicon issue on insecticidal toxins and their potential for insect pest control[END_REF][START_REF] Reynolds | Why EPA thinks resistance monitoring for Bt crops has been important and where do we go from here[END_REF]. "Behavioral resistance" to the SIT has been observed in isolated cases, where selection favors females that are able to identify and avoid mating with released sterile males [START_REF] Dyck | Sterile Insect Technique: Principles and practice in Area-Wide Integrated Pest Management[END_REF].

As with any other new control method, resistance to the lethal mechanisms of the RIDL technology or other prospective methods of autocidal biological control is therefore a potential problem. None have so far been detected in any species at various stages of developing the technology (e.g. explicitly tested in [START_REF] Gong | A dominant lethal genetic system for autocidal control of the Mediterranean fruitfly[END_REF], but that does not preclude the possibility of resistance arising after the introduction of released insects, whether through new mutations or existing mutations in some wild populations that are capable of conferring resistance.

In practice, prior to commencing any release program there is no RIDL genetic construct in the wild, so there is no positive selection pressure and any resistance-associated allele would be expected to be at very low levels (especially if it impairs fitness and therefore carries a disadvantage that would tend to reduce its frequency) unless it also conferred some unrelated pleiotropic advantage. Evidence about resistance to other pest control methods generally supports this expectation. For example, from field data for several pests, if any resistance to insecticidal Bt crops could be detected at all, allele frequencies were of the order of 10 -3 [START_REF] Gould | Initial frequency of alleles for resistance to Bacillus thuringiensis toxins in field populations of Heliothis virescens[END_REF][START_REF] Huang | Sugarcane Borer (Lepidoptera: Crambidae) resistance to transgenic Bacillus thuringiensis maize[END_REF][START_REF] Mahon | Frequency of Alleles Conferring Resistance to the Bt Toxins Cry1Ac and Cry2Ab in Australian Populations of Helicoverpa armigera (Lepidoptera: Noctuidae)[END_REF][START_REF] Stodola | Frequency of resistance to Bacillus thuringiensis toxin Cry1Ab in southern United States corn belt population of European corn borer (Lepidoptera: Crambidae)[END_REF][START_REF] Tabashnik | Delayed resistance to transgenic cotton in pink bollworm[END_REF]. This was not surprising because the toxins were derived page 7 from a soil dwelling organism Bacillus thuringiensis and so there was some prior natural exposure. However, mutant alleles associated with resistance to malathion have been found

in pinned specimens of Australian sheep blow fly Lucilia cuprina collected before the release of organophosphate insecticides [START_REF] Hartley | Amplification of DNA from preserved specimens shows blowflies were preadapted for the rapid evolution of insecticide resistance[END_REF]. This finding, in a small sample, suggests that the frequency of alleles that could confer resistance to a new product might be much higher than expected. If such polymorphisms are present at a high frequency before insecticide use, they are probably balanced polymorphisms (maintained by balancing selection such as heterozygote advantage) and therefore do not carry a clear net fitness cost in the absence of the new insecticide as would be expected for a new mutation [START_REF] Ffrench-Constant | Which came first: insecticides or resistance?[END_REF].

One possible measure to prevent or reverse the spread of resistance is to dilute the gene pool in the wild population with susceptible alleles. In radiation-based SIT, no viable hybrids are formed so there is no introgression of alleles from released insects into the wild population. With female-lethal RIDL, the F 1 males are viable and fertile so alleles from the mass-reared strain could be spread into the wild population, with the female-lethal aspect still contributing to population suppression. Alleles from bisex-lethal RIDL strains might be pushed into the wild population if there were resistance to the RIDL genetic construct so that some progeny survive and reproduce.

In previous work, we explored insect resistance to crop plants engineered to produce insecticidal toxins (Bt crops), under arrangements where refuges of non-Bt plants are provided near insecticidal crops to provide a source of alleles conferring susceptibility to page 8 the toxins [START_REF] Alphey | Combining pest control and resistance management: synergy of engineered insects with Bt crops[END_REF]Alphey et al., 2007b). We proposed that control programs could combine Bt plants with the release of Bt-susceptible homozygote insects, either natural or engineered insects, to manage economically important agricultural pests. By providing an alternative source of susceptible alleles, release of homozygous susceptible female-lethal RIDL males (or wild type insects) is predicted to slow the spread of resistance to Bt and sufficient releases can prevent or reverse the spread of the resistant allele. Here we investigate the potential for introgression of alleles from released insects to act in a similar manner, to help manage biochemically-based genetic resistance to the autocidal genetic construct itself. This could apply to any form of biochemical resistance to the lethal effector molecule expressed by the transgenic construct, for example, through alteration of the target site or increased metabolism to clear the toxin.

The evolution of allele frequencies depends on selection pressure due to relative fitnesses of the different genotypes and properties of the resistance trait. Resistance to RIDL constructs has not yet been observed, so there is no data about its potential properties, and consequently we explore a wide range of theoretical values. There is limited data about the fitness penalties of the RIDL construct in the absence of resistance. Lethality achieved by a single copy of the construct in bisex-lethal or female-lethal strains can be 97-100% (Fu et al., 2007;[START_REF] Gong | A dominant lethal genetic system for autocidal control of the Mediterranean fruitfly[END_REF][START_REF] Phuc | Late-acting dominant lethal genetic systems and mosquito control[END_REF]. Some female-lethal strains have only a modest fitness penalty, if any, in males at 0-2% lethality [START_REF] Fu | Female-specific insect lethality engineered using alternative splicing[END_REF]. Generally, from any transgenic insertion lines, only those with sufficiently close to 100% lethality in heterozygous form, and female-lethal strains with little or no fitness cost to heterozygous males, are chosen for further development. In such lines, lethality in homozygotes is 100%

or very close [START_REF] Fu | Female-specific insect lethality engineered using alternative splicing[END_REF]Stainton et al., unpublished data).

To some extent we can investigate the evolution of resistant allele frequency in isolation.

However, to assess its impact on the main objective, which is to suppress or eliminate a population, we must explore the interaction between allele frequencies and population size, and to do this we must combine a population genetic framework with a population dynamic approach. For pests or disease vectors where both sexes are responsible for damage, the number of adults emerging at each generation is a proxy measure for the level of damage.

For example, it is indicative of the level of crop damage by phytophagous agricultural pests, because every individual that survives to reproductive maturity has completed larval development on plants. For those disease vectors, such as mosquitoes, where only adult females transmit disease (when they take a blood meal), the number of adult females emerging at each generation is the relevant measure of potential harm. This might also be a good measure for some agricultural pests, such as tephritid fruit flies where oviposition damage is economically harmful.

Here, through mathematical modeling, we explore the potential for proposed release strategies to contribute to managing a hypothetical heritable biochemical resistance to RIDL genetic constructs. We examine the evolution of RIDL-resistant allele frequency and of the frequency of the RIDL construct, considering both dynamics over time and equilibria. We investigate the impact of the release ratio, initial resistant allele frequency and genetic properties of the resistant allele. First we consider dominant or recessive, page 10 completely effective resistance with no fitness costs and then gradually relax assumptions to allow the strength and fitness costs of resistance to vary across the whole range of possible values. We then explore the interaction between pest control and resistance management by modeling relative population size over time.

Methods

Parameter definitions and symbols are set out in Table 1.

Population genetics model

We use a deterministic, non-spatial, discrete-generation population genetic model of a closed homogeneous population, with random mating, and no immigration, emigration or mutation. We assume a 1:1 sex ratio.

The number of individuals surviving to reproductive maturity depends on the relative fitness values of the immature life stages. Embryonic, larval or pupal susceptibility to the lethal effect of a RIDL genetic construct is assumed to be controlled by a single autosomal locus with two alternative alleles: resistant R (frequency p in adults emerging in the wild) and susceptible S (frequency q, p + q = 1). There are three possible genotypes at this locus:

SS, SR and RR. For convenience, we use the terms "larvae" and "larval stage" to refer to all immature stages, although engineered lethality might in principle be designed to prevent egg hatch, pupation or adult eclosion instead of being fatal during larval development.

The resistant allele may have fitness costs, which we assume also take effect during the larval stage. ).

Engineered males for release ("RIDL males") carry a repressible dominant lethal genetic construct, for which they are homozygous. The RIDL construct is either bisex-lethal "BL"

(where individuals of either sex inheriting one or more copies of the construct are nonviable in the absence of the repressor), or female-lethal "FL" (which similarly renders females non-viable but allows males to survive). There are two possible alleles at the insertion site of the RIDL construct: L, representing insertion of the construct (frequency l in adults emerging in the wild), and the wild type allele w (frequency l 1 ), representing absence of the construct, giving three possible genotypes at this locus: ww, Lw and LL. We assume there is no linkage between the RIDL insertion and the locus controlling resistance and neither of these is sex-linked.

The RIDL construct imposes a relative fitness cost k H during the larval stage, where k denotes either a "target" or "non-target" of the engineered lethality. The fitness penalty H target is equal to 1 if the construct is fully lethal, and applies to both sexes for a bisex-lethal (BL) RIDL construct and to females only for a female-lethal (FL) construct. The sexspecific fitness cost H non-target imposed by an FL construct on males, either takes a lower value, to allow for a modest associated fitness penalty (for example, due to insertional mutagenesis and/or expression of a marker gene incorporated into the construct (Marrelli et page 12 al., 2006)), or is zero (having no effect). We assume that the effect of the construct is multiplicative, so the relative fitness of genotype j (j: ww, Lw or LL) in the absence of resistance is

j k j K H : 1 (1)
where j K is the number of copies of the RIDL construct that the individual possesses (0 for ww, 1 for Lw or 2 for LL).

We

J J

). Resistance may be incomplete

( 0 t RR J
), meaning that resistance might only reduce the fitness penalty associated with the L construct ( 0 ! RR J ), rather than negate its effect altogether (complete resistance, 0

RR J ).
We assume that resistance is not over-or under-dominant (

1 d d SS SR RR J J J
).

Combining all these features, the relative fitness of larvae of genotype i, j, k (i: SS, SR or RR, j: ww, Lw or LL, k: target or non-target) is then

j i k i ijk K J H \ : 1 (2)
There are potentially 18 genotypes in this model: combinations of 3 variants at the S/R locus, 3 variants at the L/w locus, and (for FL only) 2 variants in target or non-target sex.

No release. Our baseline case is where the R allele is present but not at fixation

( 1 0 0 p
) in a naïve population (i.e., one with no L alleles in the wild, 0 0 l , which would be the case before any RIDL release program commenced) and no release strategy is being employed ("NR" no release). In this scenario, only the S/R genotype is relevant (no L allele is ever present as none are introduced) and relative fitness is

i ijk \ :
for all genotypes that can arise. Using standard population genetic theory [START_REF] Gillespie | Population genetics: a concise guide[END_REF] shows that the change in resistance allele frequency in a generation, p ' , is

RR SR SS SS SR SR RR RR SR SS SR RR p pq q q p pq p p pq q pq p p \ \ \ \ \ \ \ 2 2 2 2 2 2 ) ( ) ( 2 : : : : : ' (3) Where 1 , 0 , z q p
(resistance allele present but not at fixation), the sign of p ' is determined by the sign of ) ( ) (

SS SR SR RR q p \ \ \ \
, which is never positive in our model. In the NR scenario, if there is no cost associated with the R allele it experiences no selection pressure, otherwise it would decline over time due to the effect of fitness costs.

In practice, rather than determining the R allele frequency by calculating p ' from Eq. (3), we use a simulation approach that involves calculating the frequencies of every genotype and applying the genotype-specific relative fitness values (Eq. 2). This approach allows us to analyze the model in the same way for no release (NR) and the more complicated cases of release of engineered insects (BL or FL); the latter cases involve a complex system of up to 18 difference equations that cannot be readily reduced to expressions in terms of the allele frequencies. page 14

Release of engineered insects. To model the effects of the release of RIDL males carrying either a bisex-lethal (BL) or female-lethal (FL) construct, we assume that at each generation, as adults emerge prior to mating, mass-reared adult RIDL males are released at a fixed ratio d (the "RIDL release ratio") to the number of males in the wild population (of any genotype) that survived to reproductive maturity. All engineered released males are assumed to be homozygous susceptible to the RIDL construct (SS), and homozygous for the RIDL construct (LL). Wild and released adults mix homogeneously and females mate with wild or released males at random.

Inherent in our model is a proportional release policy, where the number of insects released is proportional to the number of adult males in the current population. An obvious alternative would be a constant release policy, where the same fixed number of RIDL males is released at every generation (which could be described by the "release ratio" of that number to the number of adult males at the start of the program). Implementing the proportional policy would require more sophisticated monitoring to provide data for reasonably accurate prediction of the adult male population at each generation so that the numbers for release can be estimated. We model the proportional policy because working in terms of proportions permits us to study the evolution of allele frequencies independently of population size (so those results are independent of any assumptions about population dynamics).

Population dynamics

Density-independent populations. R 0 is the average number of progeny produced per adult pest insect that survive to adulthood in a density-independent population, commonly expressed as female offspring per adult female. For agricultural pests under field conditions, after taking into account factors such as egg hatch rates, predation and pathogens that are common to all genotypes throughout the habitat, R 0 might typically range from 3 to 11 but could be much higher in some populations [START_REF] Dyck | Sterile Insect Technique: Principles and practice in Area-Wide Integrated Pest Management[END_REF][START_REF] Dye | Models for the population dynamics of the yellow fever mosquito, Aedes aegypti[END_REF]. For NR, BL and FL strategies the (dimensionless) relative population size is determined by

t t t F R N V 0 1 2 (4) ) ( 0 1 F t t t F R F V (5)
where N t is the population size of mature adults at generation t relative to the initial pest density (set N 0 = 1), of which F t are female. The population growth rate is 2R 0 (each adult female contributes R 0 male and R 0 female larvae to the next generation) and t V and

) ( F t
V are the proportions of all offspring or female offspring, respectively, that survive to maturity (simulated by applying the fitness costs of R alleles and the RIDL construct, and the lethal effect of the RIDL construct tempered by the resistance trait, as appropriate to each genotype).

Populations subject to density-dependent processes. Density-dependent processes play an important role in the dynamical behavior of biological populations, and a number of difference equation models to portray them have been described [START_REF] Bellows | The descriptive properties of some models for density dependence[END_REF]. We adapt a general model of the form Modeling resistance to genetic control of insects page 16

t N t t e N R N 0 1 ~D (6)
where Ñ denotes absolute (rather than relative) population size, and α is a measure of the strength of density dependence. Here D

1 is related to the carrying capacity of the habitat, such that if 0

D

we recover the density-independent model, and in the limit as

f o D , 0 ~1 o t N
(the carrying capacity approaches zero and nothing survives). This system has equilibria at 0 ~* N and

D ) log( ~0 * R N . ( 7
)
To make this model dimensionless, we assume that the population starts at the equilibrium in Eq. ( 7)

( * 0 Ñ N
) and measure the population size relative to that value

1 • © § * Ñ N N t t , so that 1 * 0 N N
. It is important for our purpose to distinguish between males and females, so we express population growth in terms of progeny per female as we did in the density-independent model. Finally, we multiply by the survival proportions (V ), which implies that the genetically-induced lethality occurs after the density-dependent processes (representing a late-acting rather than early-acting genetic construct).

With parameters and symbols given in Table 1, the relative population size in our densitydependent model is determined by

t N R t t t e F R N V ) ) (log( 0 1 0 2 (8) ) ( ) ) (log( 0 1 0 F t N R t t t e F R F V (9)
To ensure that the conclusions we draw with this simple form of density dependence are valid, we have repeated our analyses using a two-parameter density-dependent function (see supplementary material).

Simulations

We assume a naïve population ( 0 0 l ) with genotypes initially in Hardy-Weinberg equilibrium (frequencies SSww: 2 0 q , SRww:

0 0 2 q p , RRww: 2 0 p ).
In the case of RIDL releases, other genotypes (Lw and LL types) arise from matings in the field.

For simplicity, we assume the RIDL construct is fully dominant and completely lethal ( H target =1) and the FL version is lethal to females but not to males; in these simulations, H non-target was set to 0. As no resistance has been detected, there are no data about its potential properties, so we consider values across the permissible range (the constraints imposed in our model are summarized in Table 1).

We previously showed that lower release ratios than typically employed for the purpose of population suppression could be useful in some cases where the main objective is to manage resistance to an independent control measure such as Bt crops [START_REF] Alphey | Combining pest control and resistance management: synergy of engineered insects with Bt crops[END_REF]Alphey et al., 2007b). However, in the context of resistance to a RIDL construct, a higher release ratio of RIDL insects would provide both more S alleles to help dilute resistance and greater selection pressure in favor of resistance, so it is not clear what the effect of a higher or lower RIDL release ratio might be.

We explore the effects of release of RIDL males by simulating different combinations of γ and ψ representing varying strength and fitness costs of resistance, over a range of release ratios from 1:10 (one RIDL male per ten males in the wild, 0.1) to 1000:1, or a subset thereof. Generally we model a release ratio of 10:1, which could be a realistic ratio for a control program; overflooding ratios quoted in practice (typically 10:1 to 40:1 [START_REF] Dyck | Sterile Insect Technique: Principles and practice in Area-Wide Integrated Pest Management[END_REF]) may be higher to allow for poor dispersal and other losses during distribution, whereas our release ratio is net of any such effects and all the released insects are assumed to be available for mating. We simulate the frequencies of all genotypes over at least 30 generations to explore transient dynamics and more as necessary (30 up to several thousands) to estimate the value of an equilibrium allele frequency. Those higher numbers are far beyond the realistic practical duration of any release program, but we use them for accurate estimation of equilibria by generating results past any transient behavior to a time when the allele frequencies should have converged.

We generally model initial R allele frequency p 0 of 0.01 to investigate resistance that is present at perceptible levels in the wild population, or 0.001 to represent a rarer allele that would probably not be readily detectable. We also use an intermediate frequency 0.1 that should be easy to detect in a field population and would cause practical concern, and some much higher values (for example, 0.7 or 0.99) to explore the potential reversal of very common resistance. A recent mutation might only have frequency of the order of 10 -6 ; we do not present results for such very rare alleles as we expect that any trend indicated by results for p 0 of 0.01 or 0.001 (and other values where simulated) would extend to lower page 19

values, and because we anticipate that it would take relatively many generations for the impact of such a very rare resistant allele to manifest.

Results

Managing resistance to a fully-lethal genetic construct

We begin by considering the baseline case where there is no resistance and then investigate the simple cases where resistance is either dominant or recessive, is complete (entirely removes susceptibility to the RIDL construct) and has no fitness costs. Later we progressively relax these assumptions to consider the full range of resistance properties and fitness costs.

Resistance is absent or lethal. There are three situations in terms of model parameters that are equivalent in that they represent the case where there is no resistance in the population: no R allele ( 0 0 p ), a candidate "R allele" that does not actually confer any resistance property (susceptibility 1

RR SR SS J J J
) and so is not a resistant R allele, or an R allele that is dominant fully lethal (relative fitness 0

RR SR \ \
) and so never confers any resistance trait. In this baseline case we consider the implications for R and L allele frequencies.

In these circumstances, a resistant R allele ( 1 RR J ) tends toward (or remains at) zero frequency, i.e. the only stable equilibrium is extinction: , the R allele's advantages will not be sufficient to achieve this and the R and L allele frequencies will settle at those minimum equilibrium values (Eqs. 10-12).

0 * p . ( 10
)
For convenience, we will refer to these minimum equilibrium values (Eqs. 10-12) as

"floor" values and to the combinations of susceptibility parameter values where they are the outcome as "floor regions".

With RIDL releases (BL or FL), a proportion ) 1 ( d d of progeny embryos (prior to the effects of R and L) will be sired by released SSLL males, and so will have genotype of the form S_L_ (meaning at least one S allele, i.e. the S/R type is SS or SR, and at least one L allele, i.e. the L/w type is Lw or LL). As release ratios would usually be 1:1 or higher

( 1 t d , with realistic 1 !! d
for a suppression or elimination program), the majority of page 22 embryos will have genotype of that form. We therefore expect that the fitness costs for SR heterozygotes (fitness SR \ ) and the effectiveness of heterozygote resistance (susceptibility SR J ) will have more influence on the evolution of R allele frequency, and hence on L allele frequency, than those of RR homozygotes do. This effect is more pronounced as the release ratio, d, increases.

Below, we confirm and expand on these general deductions by progressively relaxing assumptions about the genetic properties of resistance.

Complete dominant no-cost resistance. The strongest combination of resistance

properties is an R allele with no fitness costs (fitness 1

RR SR \ \
), conferring resistance that is complete (i.e. negates any effects of the L allele) and dominant (susceptibility 0 RR SR

J J

). Such an allele has the greatest potential to spread among the population in the presence of the RIDL construct, and simulations show that it spreads rapidly when RIDL releases are in progress (Fig. 1a,b, BL and FL). The R allele equilibrium values are above ½ for BL and above ⅓ for FL releases, which is abundant enough to allow the L allele to spread very quickly to fixation (Fig. 1c,d). FL is better than BL at controlling resistance (Fig. 1a,b) and the L allele spreads marginally faster with FL than BL (Fig. 1c,d), because the female-specific L allele is inherited through all males not just resistant males. In the absence of the RIDL construct, the R allele encounters no selection pressure and will remain constant (Fig. 1a,b, NR). page 23 Similar patterns are observed across a range of values for the initial R allele frequency.

Starting with p 0 at 0.1 (Fig. 1a,c), 0.01 or 0.001 (Fig. 1b,d The values of R allele equilibria under BL and FL strategies depend on the release ratio.

The dominant R allele has a significant advantage in the presence of the L allele and spreads as far as it can subject to dilution (which effect is stronger at higher release ratios).

As release ratio d increases, * p decreases toward ½ for BL or ⅓ for FL, and this equilibrium is approached more quickly (Fig. 2a,b). A higher release ratio also increases the rate of spread of the L allele to fixation (Fig. 2c,d

\ \

) and gives complete resistance ( 0

RR J ) but is recessive (susceptibility 1 SS SR J J
) rather than dominant, this is weaker resistance than the dominant equivalent.

The R allele only has an effect in homozygous form.

With BL release, all progeny of released males are SSL_ or SRL_ and so are non-viable, and there can be no introgression of S or L alleles into the wild population. The R allele remains at its original frequency because RRL_ genotypes (i.e. RRLw or RRLL) never arise to create selection pressure, and the L allele does not enter the population (Fig. 3, BL).

This holds true for all recessive resistance, whether or not it is complete ( 0 t

RR J

).

With FL release, all female progeny of released males are SSL_ or SRL_ and so are nonviable, but there is introgression of S or L alleles through the male line. With release ratio 10:1, the dilution effect of S alleles from released RIDL males is strong enough that the R allele declines to extinction (Fig. 3a), which means that the L allele settles at its floor equilibrium value (Eq. 12) and does so rapidly at this release ratio (Fig. 3b). Increasing the release ratio increases * l , with a limit of 0.5 as f o d .

The equilibria, and the finding that there is no change in R allele frequency with BL, are independent of initial R allele frequency.

J J

). The closer that susceptibility to the RIDL construct is to 0 (i.e. the more effective resistance is), the faster both R and L alleles spread (Fig. 4).

With FL releases, the R allele equilibrium is lower with less effective resistance (Fig. 4b).

Under the BL strategy there is very little variation in the R allele equilibrium (with releases at 10:1 * p remains close to 0.52, Fig. 4a). For both strategies, the additional survival of L alleles with stronger resistance can cause a significant increase in their frequency (Fig. 4c,d).

When we allow susceptibility to the RIDL construct to range over all permitted values

( 1 0 0 d d d SR RR J J
) with no assumptions about dominance, we see that this variation generally has more effect on * l than on * p (Fig. 5). , which is 0.01 in Fig. 5a), at realistic release ratios equilibria are independent of initial R allele frequency (we modeled p 0 of 0.01 (Fig. 5), 0.1, 0.7 and 0.99 with d = 10), and where resistance is recessive ( 1 SR J

) the R and L alleles tend to their floor values (Eqs. 10-12) (Fig. 5b-d). Where resistance is complete and dominant (i.e., one or more copies of the R allele prevent any susceptibility to the genetic lethality 0 RR SR

J J

), the R allele spreads as far as it can subject to the effects of dilution (Fig. 5a,b) and the L allele goes to fixation (Fig. 5c,d). The L allele also goes to fixation for parameter values in a region around that corner of parameter space ( SR J and RR J near zero), where the lethality of the RIDL construct is still overcome by the very effective resistance compared with the effects of dilution. What happens in between is slightly different for BL and FL strategies.

In the BL scenario, for any given RR J , as heterozygote susceptibility SR J increases toward higher values (weaker heterozygote resistance), * p increases (Fig. 5a) and * l decreases (Fig. 5c) toward extinction (its floor, Eq. 11). This is the result of fewer SRL_ individuals surviving; this reduces the effect of dilution by S alleles and introgression of L alleles from released RIDL males because that is the genotype of their potentially resistant offspring.

For given SR J , as homozygote susceptibility RR J increases, * p (Fig. 5a) and * l (Fig. 5c) both decrease (although not by much) because reduced survival of RRL_ types reduces the frequency of R alleles a little. With lower release ratio (d = 5), the equilibrium surfaces are similar, but with higher * p values due to less dilution, and lower * l values, with L allele reaching fixation in a slightly smaller region of parameter space near the corner 0 RR SR

J J

than it did with d = 10. At even lower release ratio (d = 1) those results remain robust, although the * p surface appears a little less flat near one edge ( 0

RR J ),
and again, the equilibrium results are the same with p 0 0.01 as with 0.99.

In the FL case, the outcomes are more sensitive to some parameter values. As in the BL case, for any given SR J , as homozygote susceptibility RR J increases, the equilibria * p and * l decrease a little (Fig. 5b,d), and over most of the range of values, holding RR J fixed and increasing SR J (less effective heterozygote resistance) reduces * l towards its floor (Eq. 12) (Fig. 5d). This is consistent with fewer target SRL_ females surviving and so reducing the introgression of L alleles from released RIDL males through their female progeny. There is distinct interesting behavior in one corner region, where 0 RR J (homozygotes have complete resistance and so are as fit as SS genotypes), and SR J is near 1 (resistance is closer to recessive than dominant); this situation where 0 RR J is investigated fully in the supplementary material. With release ratio 10:1, there is little change in the R allele equilibrium frequency * p over most of the surface (Fig. 5b). Holding RR J constant, increasing SR J can slightly increase or decrease * p at different values of SR J , and again there is atypical behavior in the same region. With a release ratio 10:1, the * p and * l equilibria surfaces are the same regardless of initial R allele frequency (identical equilibrium results for p 0 0.01 and 0.99). With a lower release ratio (d = 5 or 1, supplementary Fig. S1), the surfaces appear similar to Fig. 5b&d except that the * l floor is lower (Eq. 12) and the * p and * l floor regions extend further from the edge where 1

SR J ;
the lower release ratio, the greater the observed effects. The outcome in the atypical corner is affected by initial R allele frequency; when p 0 is higher there are more pronounced peaks in the region of interest (see supplementary material results).

The impact of fitness costs of resistance.

Where the R allele has fitness costs, recessive resistance ( 1

SR

J

) is no longer a special case for the BL strategy; the costs cause the R allele to decline to extinction rather than remain unchanged at its initial frequency. The existence of fitness costs for resistance to a FL RIDL construct removes the interesting behavior discussed above from the region of parameter space where SR J is near 1 (resistance is nearly recessive) and 0 RR J (homozygote resistance is complete). (We modeled p 0 of 0.01, 0.1 or 0.8 with d = 10, and p 0 of 0.01 with d = 1, to confirm both these findings). As a result, the BL and FL equilibrium surfaces against susceptibility to RIDL are similar in shape (but not scale) and behavior, so we illustrate only FL. Resistance with high heterozygote susceptibility SR J is not sufficiently effective, in the presence of fitness costs, for the R and L alleles to persist above their floors (Eqs. 10-12)

and those floors extend over a larger region (mostly defined by a SR J threshold).

A resistant allele that is very costly (high fitness penalty in heterozygotes and homozygotes) would be harder to control when it is very abundant than when it is rare.

However, such a high-cost allele would not be expected to be abundant in practice, because of the selection pressure. With 10:1 release ratio, p 0 only makes a difference to * p and * l for the most costly R alleles ( SR \ and RR \ near zero); for such alleles, higher p 0 gives a smaller floor region (the edge occurs at higher SR J ) but does not alter equilibrium values elsewhere. That is, for costly resistance with given susceptibility, p 0 determines whether the R and L alleles will persist above their floors (Eqs. 10-12) but, if they do, it does not change the equilibrium values at which they will persist. These results suggest some kind of bifurcation of surfaces (floor versus persistent equilibria) for each of the R and L alleles for a given release ratio, with p 0 determining on which surface the result lies. The critical p 0 appears to be moderated by the release ratio and fitness cost of resistance, and if it does not fall between 0 and 1, floor values will result either at all values of susceptibility or only for recessive resistance (the edge where 1

SR J

). At the edge of the floor, the transition can be continuous (Fig. 6a,c) or discontinuous (Fig. 6b,d).

A lower release ratio (1:1 instead of the 10:1 shown in Fig. 6) reduces the SR J threshold thereby enlarging the floor regions (for example, with d = 1 and costly resistance, floors for both R and L alleles extend across the whole range of susceptibility values for BL and FL).

This occurs because there is less selection pressure in favor of the R allele at lower release ratios. With effective resistance (on the low susceptibility SR J side of that threshold) the L allele goes to fixation over a smaller range of susceptibility values and * p is higher. The R allele persists at higher equilibrium value because the lower release ratio lessens the effect of dilution by S alleles from released males.

These effects of reducing the release ratio suggest that for any given set of parameter values (γ, ψ) and initial R allele frequency there is a critical release ratio, below which the R allele tends to extinction and the L allele tends to extinction (BL) or its non-zero floor (FL). This is borne out by simulations. This critical release ratio depends on the genotypes' relative fitness values, with regard to both costs and efficacy of resistance and, for some parameter combinations, the initial R allele frequency. For example, with 9 . 0 . Even so, as results in the next section show, the spread of a page 32 resistant allele will not necessarily have a significant impact on suppression of the population.

Population control and resistance management

For practical purposes, it is not the frequency of a resistant allele that is directly of interest, but the size of the pest population. We explored the interaction between the evolution of allele frequencies and the impact on the relative size of the pest or vector population. We also examined the relationship with the insect's net reproductive rate R 0 (the per capita rate of population growth in a generation), employing a higher release ratio for a pest with higher growth potential. We limited these investigations to situations where the R allele will not decline to extinction, because that is where concern would focus, and compared the outcomes with the equivalent release regimes in the absence of any resistance. We selected release ratios that would lead to a successful program (local elimination of the pest) in the absence of resistance.

With no release and bisex-lethal RIDL release strategies, the females are always half of the population. Under the female-lethal RIDL release strategy this is not true; the relative number of females is less than half the relative total population size because females are targeted and are killed off disproportionately. The effect of resistance is qualitatively similar whether measuring the total number of insects or the number of females (illustrated in supplementary Fig. S3, compare panels a & d).

Under our simulated proportional release policy, female-lethal releases appear to perform better in terms of relative population size than bisex-lethal releases at the same release ratio (Figs. 7,S3 & S4). This occurs as the number of insects released is proportional to the number of males emerging in the wild; the BL strategy kills males and females in equal numbers but FL strategy kills only females. The same release ratio for both strategies involves more insects being released under the FL strategy than the BL strategy, so the two are not directly comparable.

Density-independent populations. The presence of a resistant allele that spreads to nonzero equilibrium has greater effect on the population size if it is a stronger R allele (with low fitness costs and low susceptibility to the lethal genetic construct) and/or if the insect has a high per capita population growth rate. A weaker R allele (with a smaller reduction in susceptibility and higher fitness costs) can have a fairly small impact on population control, by marginally slowing the rate at which releases drive the pest to local elimination (supplementary Fig. S3a). If the pest has a low population growth rate, such resistance can have negligible impact on population management (with lower R 0 and d values, as used in Fig. S3b, the results with and without resistance are indistinguishable on the equivalent of Fig. S3a until the relative population size is below 10 -20 ). A stronger R allele has a greater effect on population size, although with lower R 0 the population may remain suppressed below its original size for many generations (Fig. S3b). At worst, at a release ratio that would be very effective in the absence of resistance, strong resistance in a pest with higher R 0 can mean that instead the releases only slow the rate of explosion of the population (Fig. S3c). In all instances we modeled, releases did at least as well as no release: even complete dominant no-cost resistance (fitness 1

RR SR \ \ , susceptibility 0 RR SR J J
) resulted in significantly slower population growth than for the uncontrolled population.

As described above, the initial R allele frequency does not make a difference to the equilibrium R or L allele frequencies for most combinations of fitness costs and susceptibility. If the allele starts at a much higher frequency its impact on the population size in the early generations could potentially affect the relative population size after the equilibria have been reached. However, we generally observed only a slight slowing of population decline or increase of population explosion (by a small number of generations).

This might or might not be large enough to have a significant economic or public health effect. For example, with parameters as in Fig. S3, the largest difference was with the stronger resistance with the higher growth rate (Fig. S3c), where the population became more than a million times its original size after 14 generations from 9 . 0 0 p , compared with 19 from 01 . 0 0 p . The impact was much smaller for the other panels.

Populations subject to density-dependent processes. Broadly, we draw the same general conclusions about the impact of resistance to the RIDL construct when there are densitydependent processes operating as for density-independent populations. The main difference is that the population is assumed to start at equilibrium, so the relative population size will remain at 1 with no releases (Fig. 7). The RIDL construct fitness penalty takes effect after density-dependent mortality, because of its late-acting lethality, and successful RIDL releases can drive the population to local elimination. Where resistance disables the efficacy of the RIDL construct, the population will return to equilibrium densities.

Density-dependent processes do not alter the prediction that a spreading resistant allele has greater effect on the population size if it is a stronger (low susceptibility and low cost) R allele and/or if the insect has a high per capita population growth rate. The presence of a weak R allele can mean releases take slightly longer to drive the pest to local elimination (Fig. 7a); with low population growth rate the impact on population control might be negligible (e.g. as Fig. 7a except with the lower R 0 and d values). A strong R allele is more detrimental for the effectiveness of releases, although with lower R 0 the initial suppression, while resistance is still uncommon, is enough to keep the population below its original size for many generations (Fig. 7b). At high release ratios, strong resistance spreading in a pest with higher R 0 might only suppress the pest for a few generations before the population recovers to near natural equilibrium densities (Fig. 7c). If that R allele is initially widespread (e.g. as Fig. 7c but with 9 . 0 0 p instead of 0.01) that initial suppression phase does not occur and pest densities move quickly to that new level just below the natural equilibrium.

We obtained qualitatively similar results with a two-parameter density-dependence process (see results in supplementary material and Fig. S4). A weaker R allele (with high susceptibility to the genetic lethality and some fitness costs) can make negligible difference to the effectiveness of controlling a low growth rate population, and the impact of resistance increases with the strength of the R allele and the population growth rate.

In all instances modeled, we found that releases reduced the population below the natural equilibrium. Even complete dominant no-cost resistance (fitness 1

RR SR \ \ , susceptibility 0 RR SR J J
) resulted in population densities slightly below equilibrium and could exhibit suppression temporarily while resistance was uncommon. (For example, with other parameters as in Fig. 7b i.e. p 0 0.01, d 10 and 2R 0 6, population size declined to 0.03 relative to equilibrium in 3 (BL) or 4 (FL) generations, then recovered to 0.45 (BL) or 0.41 (FL) by the 17th generation.) The temporary period of effective suppression can last longer if starting with a rarer allele.

Discussion

We have shown that a RIDL approach could potentially limit the effects of biochemically based resistance against the RIDL genetic construct (L allele) engineered into the released strains, due to a built-in dilution mechanism. Resistant alleles are unable to spread to fixation, although they can become more common than the alternative susceptible allele.

However, the evolution of resistant allele frequency is subject to a complicated mix of competing forces.

A resistant allele with no fitness costs has a net advantage in a population containing the RIDL genetic construct, the magnitude of which depends on the effectiveness and dominance of resistance (which we express in terms of genotype susceptibility to the RIDL construct). That advantage depends on the frequency of the RIDL allele, which in turn depends on the release ratio. However, whilst a higher release ratio increases the frequency of the L allele and thereby acts in favor of resistance, it also increases the dilution of the gene pool by susceptible alleles to reduce resistance. The existence of fitness costs associated with resistance adds another complication by affecting the direction and page 37 magnitude of the resistant allele's net advantage or disadvantage. This is a multidimensional problem for which several parameter values cannot reasonably be predicted because resistant alleles are likely to be rare or non-existent at the commencement of a release program, and may be difficult to measure or estimate accurately even if resistance has been identified, especially if resistance is recessive or nearly recessive so that it is only observed in homozygous form.

In the absence of resistance, increasing the ratio of released RIDL males to males in the wild has a beneficial impact on the size of the pest population. As with classical SIT, there is a critical release ratio above which the size of the population will be reduced [START_REF] Dyck | Sterile Insect Technique: Principles and practice in Area-Wide Integrated Pest Management[END_REF][START_REF] Phuc | Late-acting dominant lethal genetic systems and mosquito control[END_REF][START_REF] Thomas | Insect population control using a dominant, repressible, lethal genetic system[END_REF]. Here we have identified that there is another critical release ratio below which resistance will be prevented from spreading through the population and the RIDL construct will not spread among individuals of the target sex(es).

Clearly, there may be situations where those two requirements are incompatible, which arise when resistance is more effective and/or has lower fitness costs. However, the population might still be effectively controlled even if resistance becomes common, and with weak resistance that is initially rare there may even be little practical impact for population control. Stronger, fitter resistance has greater potential to be detrimental for the aim of reducing population size. Some of the low release ratios that we have modeled would not intentionally be used in an actual population control program. However, these values might be relevant outside the program's target area if there are regions to which engineered insects could emigrate. It is known that, with density-dependent populations, SIT releases can in theory have unintentional consequences -sometimes detrimental -in neighboring areas, although a related modeling study suggests this would not occur with (bisex-lethal) late-acting RIDL strategies [START_REF] Yakob | Aedes aegypti control: the concomitant role of competition, space and transgenic technologies[END_REF]. Our findings here suggest that any resistance to the genetic construct in areas of low actual release ratio would be less likely to spread, because of the page 39 lower favorable selection pressure; it might therefore be expected to interfere less with population reduction outside the target region, than inside it. This is unlike resistance to chemical insecticides, where lower doses (such as spray drift, or early prototype "mediumdose" GM crops) tend to encourage resistance by allowing survival of partially resistant individuals. A model incorporating spatial aspects would be required to investigate this matter fully.

We have not simulated the situation where the RIDL genetic construct is not fully lethal ( H target < 1). We might reasonably expect a trade-off between population suppression and resistance management, but the competing selective forces acting on the evolution of resistance might lead to counterintuitive results, so this could be worthy of exploration. We did not simulate scenarios where the female-lethal RIDL construct has a modest fitness cost in males ( H non-target > 0). This also could be potentially interesting in conjunction with the population dynamics.

We assumed that the insects released are homozygous susceptible. The dynamics could be different if quality control processes are not able to ensure this. In particular, if resistance is of low dominance, it might be that the release population includes a small number of heterozygous susceptible RIDL males (SRLL types as well as SSLL). In principle our model could be adapted to explore this issue.

Compared with the proportional release policy of our model, a constant release policy does not enable population genetics to be studied separately from the population dynamics, so page 40 some population growth assumptions have to be made, however simple. If the release ratio of released males to initial population of males, say, is high enough to reduce the population initially, in subsequent generations the ratio of released males to actual males in the wild increases. This could mean that the resistance becomes more able to spread as the effective release ratio increases, which in turn could impact the population size. This adds a further layer of competing influences above those encountered in our proportional policy model. In practice, what is implemented might be somewhere between those two, with a larger release quantity at first, reduced in one or more steps as the population declines.

Some mosquito populations vary with environmental factors such as rainfall. Seasonal variation could be modeled by varying the carrying capacity over time. With our proportional release policy, the population genetics are independent of the population dynamics, so our conclusions about allele frequency evolution would not change (unless the resistant phenotype was environmentally dependent, for example, with reduced ability of eggs to survive desiccation). With a constant release policy, the effective release ratio would vary (being at its minimum when the population is at its seasonal peak, and conversely), and seasonality would vary the balance of competing selective forces.

Releasing fertile, susceptible, wild type insects would allow introgression of susceptible alleles into the population in a similar manner to release of RIDL insects. However, given the difficulty of separating wild males from females on an industrial scale, both sexes would probably have to be released and any released fertile females would contribute to population growth. Susceptible alleles would be inherited through both sexes and none of page 41 the progeny would suffer reduced survival (unlike progeny of RIDL males), so this should be more effective at reducing the frequency of resistant alleles than equivalent RIDL releases. We could incorporate wild type releases into our models, to explore combinations of release strategies to manage resistance and maintain the efficacy of an established RIDL release program. In theory, we could also consider the sequential or alternating use of bisex-lethal RIDL releases and female-lethal RIDL releases, with one being used to combat unrelated resistance to the other while still managing the population size. This would only be appropriate if the two genetic constructs involved had independent modes of action and there was no cross-resistance. Given the molecular basis of the RIDL strains currently in existence or under development [START_REF] Fu | Female-specific insect lethality engineered using alternative splicing[END_REF][START_REF] Fu | Female-specific flightless phenotype for mosquito control[END_REF][START_REF] Gong | A dominant lethal genetic system for autocidal control of the Mediterranean fruitfly[END_REF][START_REF] Phuc | Late-acting dominant lethal genetic systems and mosquito control[END_REF]Stainton et al., unpublished data;[START_REF] Thomas | Insect population control using a dominant, repressible, lethal genetic system[END_REF], these may or may not be reasonable assumptions depending on the particular combination of BL and FL strains used.

Our model could be modified to investigate further questions relevant to resistance management strategies. Any resistant allele would likely be very rare in a natural population, and stochastic effects would be expected to have implications at low frequencies. For example, with a stochastic model we could explore the probability that a resistant allele will go extinct within a specified number of generations. "Behavioral resistance" (see Introduction), rather than biochemical resistance, would require us to replace our assumption of random mating, which would potentially result in different conclusions. We note that a broadly similar model has also been used to simulate the intentional spread through a population of an allele resistant to a lethal gene, as a means of page 42

driving an anti-pathogen gene through a vector population, with a "killer" gene conferring an advantage to a "rescue" gene to which would be linked a genetic construct that inhibits transmission of a target pathogen [START_REF] Gould | A Killer-Rescue system for selflimiting gene drive of anti-pathogen constructs[END_REF].

A potential concern about releasing GM insects into the wild is that the inserted DNA may have unforeseen consequences. This has particularly been raised regarding strategies that aim to replace a wild population, for example with a version engineered to be unable to transmit a pathogen. One advantage claimed for genetic SIT strategies over population replacement strategies is that autocidal engineered insects are programmed to die and therefore the lethal genetic construct should die out if releases cease. This relies on the construct having some fitness costs even if it does not retain its efficacy. We have shown that resistance to the construct might significantly increase its frequency. In theory, there are combinations of genetic properties of resistance for which the construct could eventually be driven to fixation; that scenario is highly unlikely to be played out in real life because the release program could be stopped if substantial resistance were detected, long before the RIDL construct became common, and the wild type would reinvade through immigration.

It is remarkable that a fairly simple two-locus genetic system with forcing (from releases) can produce such a complicated array of potential behaviors. Some combinations of genetic properties of resistance can generate unexpected results, especially with the femalelethal strategy, and small changes in a single parameter can make the difference between resistance going extinct and resistance spreading so widely as eventually to render the lethal construct largely ineffective. However, the theoretical equilibrium results are not always relevant to planning suppression or elimination programs, where it is short to medium term behavior (the transient dynamics) that is important. Ultimately, it is the effect of the resistance on the population size that matters and that can vary from negligible to substantial impact on the results of a release program, with strong (effective and low-cost)

resistance and high growth-rate insects being harder to control. Fortunately, those stronger resistance alleles are the ones that are easier to detect in the population either during preparatory trials or during early stages of a release program.

Under an envisioned RIDL release program, it would be important to detect any resistance in the pest population that was capable of having a significant detrimental impact on the program's effectiveness in time to take remedial action. Our population dynamic models predict a clear change in the pattern of pest densities (for example, Figs. 7c & S3c) as such resistance begins to spread. A properly designed monitoring strategy would be alert to such changes. Wild-caught individuals could then be used to screen for cross-resistance against a panel of available RIDL strains that use different lethal mechanisms, so that operations (or further trials, if this occurred before full implementation) could switch to deploying an appropriate alternative RIDL line before significant lasting loss of efficacy could result. page 50

J J

). The L construct is fully-lethal to targets, with no other fitness costs:

H target = 1 and H non-target = 0. The release ratio is 10:1 ( 10 d

). Initial R allele frequency

Tables

  ), all yield the same broad pattern and identical equilibrium frequencies. Even from near fixation with the R allele declining rapidly to that value and the L allele spreading very rapidly from 0 to fixation.The shapes of these allele frequency curves against time are examples of the patterns that can be seen more generally, not only for complete, dominant, no-cost resistance. The first two generations of RIDL releases (BL or FL) into a naïve population may exhibit distinctive changes (e.g. clearly visible in Fig.1dFL), because those are when Lw (first generation) and LL (second generation) individuals first appear among the population emerging in the wild. Apart from this, the frequency of a spreading R allele is consistent with the sigmoid pattern often associated with resistance generally: a potentially long period (for a very rare allele) with little change, followed by a period of rapid increase, then leveling off approaching the equilibrium frequency (Fig.1a,b, BL, FL). The R allele never reaches fixation with RIDL releases because of the dilution effect of adding new S alleles at each generation, but the equilibrium frequency can potentially be close to 1 (see later examples). The same general patterns can be observed in the frequency of a RIDL genetic construct (L allele), except that the L allele can reach fixation with sufficiently fit and effective resistance (Fig.1b,d), and if the R allele starts at high frequency the L allele can at first increase to a higher level and then decline to its equilibrium value. page 24

  that the alleles reach their equilibria more quickly, especially with the lower release ratios. The limits for * p (½ with BL, ⅓ with FL) can be deduced from the genetics (see supplementary material Appendix B). Complete recessive no-cost resistance. If the R allele has no fitness costs (fitness 1 RR SR

  Fig. 3 illustrates p 0 of 0.001. With p 0 0.1 broadly similar patterns over time are observed. Even with very high p 0 0.99 the equilibria (or no change) are the same, but the transient dynamics differ. With FL releases, R frequencies decline monotonically from that high level and the L allele frequency increases at first (reaching a maximum exceeding 0.62 in the 15th generation) before declining to the same equilibrium value. Incomplete no-cost resistance. Resistance might reduce the fitness penalty of the RIDL construct rather than cancel it, i.e. resistance might be incomplete. The effectiveness of a dominant R allele with no fitness costs (

  page 29With both BL and FL strategies, as fitness costs associated with the R allele increase, especially those for SR heterozygotes, the R allele goes extinct for a greater range of heterozygote susceptibility values, and where it persists it tends toward a lower equilibrium(compare Fig. 5b,d no costs, Fig. 6a,c low costs and Fig. 6b,d high costs). This reduction in * p is smaller for more effective R alleles (low values of susceptibility SR J and RR J ).

  recessive incomplete resistance with nearly dominant, moderately high fitness costs), the critical release ratio for FL is roughly 37effective resistance and/or lower fitness costs give a lower critical release ratio, which could be significantly below the release ratios expected for a suppression or elimination program.Homozygote costs and susceptibility have little effect compared with heterozygote costs and susceptibility. For example, having moderate or high homozygote fitness in homozygotes is unimportant and it is mainly the phenotype in heterozygotes that determines whether resistance spreads.Regardless of p 0 , some combinations of parameter values can give rise to a large discontinuity at the critical release ratio while others do not. For example, and they become so fairly quickly -p reduces from 0.9 to just above 0.05 in 30 generations) a release ratio a little over the critical value need not be cause for practical concern

  Stronger resistance (with low fitness costs and low susceptibility to the lethal genetic construct) is the main threat to the success of a genetic control release program. Such resistant alleles have greater selection pressure in their favor and tend to spread very rapidly, even from low levels. At worst, against target insects that have little or no regulation by density-dependent processes, the presence of such resistance should therefore become apparent quickly, as monitoring data would reveal a clear detrimental change in the rate of population growth after a promising start. Similarly, population monitoring should identify such resistance among populations that are subject to density-dependent mortality, as a new equilibrium density would be reached fairly rapidly. In practice, though, preimplementation preparatory trials should be expected to detect such fit and effective resistant alleles, unless they are very rare. All RIDL strains incorporate a heritable genetic marker that allows easy identification, so survival of RIDL heterozygotes (Lw individuals) should be readily detected from trapping data during release trials. Laboratory trials involving crosses with the wild population might potentially identify resistant alleles, present in the field, at a much earlier stage. The detection threshold could be of the order of 0.01 allele frequency or lower, depending how recessive the resistance trait is (a resistance effect in SR heterozygotes is easier to detect).

Figure 1

 1 Figure 1

Figure 2

 2 Figure 2

Figure 3 Figure 4

 34 Figure 3

Figure 5

 5 Figure 5

Figure 6

 6 Figure 6

  R allele has fitness costs: relative fitness (a & c

Figure 7

 7 Figure 7

  Figure1

Figure 7

 7 Figure 7

  The relative fitness of immature individuals of genotype i (i: SS, SR or RR) is

	i \ , where fitness is measured relative to SS homozygote wild type insects (so	\	SS	1	).
	Fitness costs are not over-or under-dominant (	\	RR	d	\	SR	d	\	SS	1

  page 20As we have no mutation or immigration in our model, no resistant allele ever arises if it is initially absent. A dominant fully lethal allele will not survive into subsequent generations.

	present in the wild population (	p	0 !	0	,	J	RR	1	and	\	SR	!	0	). The equilibria set out
	A candidate allele that turns out to be non-resistant ( above for no resistance (Eqs. 10-12) are minimum values for allele equilibria in the 1 RR J ) is not of interest. For
	completeness, note that such an allele has no fitness advantage, and will either remain presence of resistance; they will still apply unless the R allele's properties are such that the
	stable (under NR or BL scenarios if it has no fitness costs) or decline to extinction due to fitness advantages outweigh the fitness penalty of the L allele and the effects of dilution by
	dilution and/or fitness costs (in all other situations). S alleles (which depend on the release ratio d). This can occur if the R allele is not lethal
	(fitness RR \ and perhaps SR \ not zero) and is effective against the RIDL construct
	(susceptibility RR J and perhaps SR J below 1). If so, the R allele can persist (	p	* !	0	), and
															l	!	(NR,
	BL) or	l !	1 ½	d d								\	SR	\	RR	and
	J	SS	J	SR	J	RR	1
								BL:		0
	Appendix A) to be
								FL:	l	1 ½	d d

The evolution of L allele frequency differs between BL and FL strategies. With no resistance and BL releases carrying a fully lethal L construct, all progeny of released males are non-viable, so no Lw or LL genotypes survive and the L allele does not spread. (There is no introgression of other alleles either, including any S alleles). Frequency l remains zero in all generations (frequencies are measured in adults emerging in the wild, and so exclude the released males). * l (11) Under the FL strategy, in the absence of resistance all Lw or LL females die (and so all adult females are ww), but males survive and so the L allele will introgress into the population. The equilibrium L allele frequency can be shown (see supplementary material * (12).

This is a stable equilibrium that increases with increasing release ratio.

Resistance is present. We now consider the situation where a viable resistance trait is so some target LL or Lw genotypes can survive (both sexes are targeted by BL and females are the target for FL), which allows the L allele to persist at higher frequency ( 0 * * (FL)). In a region of parameter space near 0

  As predicted, varying SR susceptibility has a greater effect on equilibria than varying RR susceptibility; the slopes of

	the * p and * l surfaces are generally steeper for fixed RR J than for fixed SR J (this is more
	visible for * l ). Except for the special case of recessive resistance (	J	SR	1	) with BL (where
	p	*	p	0

Table 1 Parameters and symbols

 1 number of copies of the RIDL construct 0 for ww, 1 for Lw or 2

	Symbol Parameter

j K t F

relative population size of mature females at generation t (see Eqs. 5 or 9) t V simulated proportion of all offspring that survive to maturity

Acknowledgements

This work was supported by the Royal Society (through a Royal Society University Research Fellowship to MBB), the UK Biotechnology and Biological Sciences Research Council (BB/C512702/1, MBB; BBS/S/J/2005/12055, NA) and by a grant to the Regents of the University of California from the Foundation for the National Institutes of Health through the Grand Challenges in Global Health initiative (LA).

Role of the funding sources

None of the funding sources had any role in study design, the collection, analysis or interpretation of data or results, the writing of the manuscript or the decision to submit it for publication.

Disclosure statement

LA is an employee and NA was a collaborative DPhil/PhD student of Oxitec Limited, which therefore provided salary and other support for the research program. Those authors have shares or share options in Oxitec Ltd (LA) or an indirect interest therein (NA). Both Oxitec Ltd and the University of Oxford have one or more patents or patent applications related to the broader subject of this paper.

Contributors

All authors designed the study and analyzed the results. NA performed the detailed analysis and simulations and wrote the manuscript. MBB and LA edited the manuscript.

All authors have approved the final article.

Supplementary Figure S1

Supplementary Figure S4


Laura Fumanelli
email: lfumanelli@fbk.eu

Mimmo Iannelli
email: mimmo.iannelli@unitn.it

Ahmed Hussnain

Olivier Janjua
email: janjua@science.unitn.it

Jousson
email: jousson@science.unitn.it

Hussnain Ahmed Janjua

Olivier Jousson

Mathematical modeling of bacterial virulence and host-pathogen interactions in the Dictyostelium/Pseudomonas system

Keywords: Host-pathogen dynamics, Amoeba-bacteria interaction, Population models, Bi-stability

We present some studies on the mechanisms of pathogenesis based on experimental work and on its interpretation through a mathematical model.

Using a collection of clinical strains of the opportunistic human pathogen

Pseudomonas aeruginosa, we performed co-culture experiments with Dictyostelium amoebae, to investigate the two organisms' interaction, characterized by a cross action between amoeba, feeding on bacteria, and bacteria exerting their pathogenic action against amoeba. In order to classify bacteria virulence, independently of this cross interaction, we have also performed killing experiments of bacteria against the nematode C. elegans.

Introduction

Virulence is a complex and often multi-factorial phenotypic trait that quantifies the degree of pathogenicity of an organism, by establishing its ability to mount an infection and to cause a disease. Genetic determinants of bacterial virulence are potential drug targets, and represent an alternative to the generic cellular and molecular mechanisms targeted by most antibiotics used to date for the control of infections in humans and animals. Model host organisms are used to identify and study the function of virulence factors in human pathogenic bacteria. Mammalian hosts are preferentially used for this purpose, because they are thought to mimic the situation occurring in an infected human patient. However, a number of non-mammalian organisms, such as the amoeba Dictyostelium discoideum, the nematode Caenorhabditis elegans, the insect Drosophila melanogaster and the fish Danio rerio have been recently used as alternative hosts, with the advantage of avoiding ethical problems and allowing us to assess bacterial virulence in an easier way than by using mammalian hosts.

Dictyostelium amoebae are soil microorganisms that feed on bacteria and have been utilized to assess the virulence of many different bacterial genera, including Klebsiella [START_REF] Benghezal | Specific host genes required for the killing of Klebsiella bacteria by phagocytes[END_REF], Aeromonas [START_REF] Froquet | Alternative host model to evaluate Aeromonas virulence[END_REF] or Pseudomonas [START_REF] Cosson | Pseudomonas aeruginosa virulence analyzed in a Dictyostelium discoideum host system[END_REF][START_REF] Pukatzki | The human pathogen Pseudomonas aeruginosa utilizes con-served virulence pathways to infect the social amoeba Dictyostelium discoideum[END_REF]. In this model system, pathogenic bacteria defend against predation by Dictyostelium, activating virulence pathways that impede amoeba cells to grow and to form phagocytic plaques. The system has been almost exclusively used to assess the role of specific genes in virulence, by comparing the virulence phenotype of mutants versus wild-type strains. To our knowledge, the Dictyostelium model has not been applied to date to determine the natural virulence range of pathogenic bacteria obtained from clinical samples. Mathematical models concerning bacterial growth have been recently developed within different contexts motivated by nosocomial infections, though, to our knowledge, no specific model has been considered for describing the particular amoeba-bacteria interaction such as described above. In particular, some work has been devoted to modeling the interaction between pathogens and cells of the immune system [START_REF] Nowak | Virus dynamics[END_REF] or between bacteria and phagocytes, including antibiotic resistance [START_REF] Austin | The dynamics of drug action on the within-host population growth of infectious agents: melding pharmacokinetics with pathogen population dynamics[END_REF][START_REF] Boldin | Relative effects of barrier precautions and topical antibiotics on nosocomial bacterial transmission: Results of multi-compartment models[END_REF][START_REF] D'agata | Modeling antibiotic resistance in hospitals: The impact of minimizing treatment duration[END_REF][START_REF] D'agata | The Impact of Different Antibiotic Regimens on the Emergence of Antimicrobial-Resistant Bacteria[END_REF][START_REF] Imran | The dynamics of bacterial infection, innate immune response and antibiotic treatment[END_REF][START_REF] Webb | A model of antibiotic resistant bacterial epidemics in hospitals[END_REF].

In this study, we have investigated the dynamics of virulence both per-forming experimental measures and building a mathematical model for interpreting the results. We have used a collection of clinical strains of the opportunistic human pathogen Pseudomonas aeruginosa to perform co-culture experiments with Dictyostelium. The interaction between the two organisms is characterized by a cross action between amoeba, feeding on bacteria, and bacteria exerting their pathogenic action against amoeba. In order to classify bacteria virulence, independently of this cross interaction, we have also performed killing experiments of bacteria against the nematode C. elegans.

A mathematical model was developed to describe the growth of the amoeba-bacteria system, taking into account the specific features underlying the interaction. We discuss the behavior of the model with respect to a basic parameter representing the virulence of the bacterium. The theoretical results are compared with the experimental ones. The model does not fall within the class of traditional prey-predator models because not only does an amoeba feed on bacteria, but it is in turn attacked by bacteria that kill its cells, so that the model must include a feedback term modeling this further aspect of the interaction. The resulting model shows existence of multiple steady states and bi-stability that seems to be responsible for the observed behavior.

Virulence assays

The focus of this study is the investigation of the mechanisms underlying the interaction of different virulent strains of P. aeruginosa with D. discoideum NC4; nevertheless a first set of experiments with C. elegans N2 as a target of the pathogenic process has been performed, in order to get a measure of the pathogenic action of P. aeruginosa independently of the predation exerted by D. discoideum. In fact in our context the nematode C. elegans N2 can be considered completely passive because during the time of the culture it does not feed on P. aeruginosa; thus the killing percentage that we obtain in the culture can be assumed as a measure of virulence, allowing us to order the different strains of P. aeruginosa according to their pathogenic strength.

Strains and growth media

A set of 27 P. aeruginosa clinical strains, isolated from patients suffering from urinary tract infections, wound ulcer and respiratory tract infections, were used in the experiments. Pseudomonas aeruginosa strains were obtained from the microbiology departments at S. Chiara Hospital of Trento (Italy) and at the University Hospital of Verona (Italy). Growth media were LB for P. aeruginosa culture and maintenance, SM broth for D. discoideum NC4 using Klebsiella aerogenes as a food source, and NGM agar for C. elegans using E. coli OP50 as a food source.

Virulence assays with C. elegans N2

Virulence of different strains of P. aeruginosa was measured independently of its interaction with D. discoideum by performing killing assays with C. elegans.

Table 1

The C. elegans fast killing assays [START_REF] Tan | Killing of Caenorhabditis elegans by Pseudomonas aeruginosa used to model mammalian bacterial pathogenesis[END_REF] The percentage of killing ranged from 6.6% to 100% (Table 1). No mortality of C. elegans was observed using the non-pathogenic strain E. coli DH5α as a control. Among the 27 P. aeruginosa specimens tested, 15 were highly virulent (80-100% C. elegans killing), 7 strains were moderately virulent (50-80% killing) and 5 were weakly virulent or non-virulent (less than 50% killing).

Virulence assays with D. discoideum NC4

Measurements of D. discoideum growth on a lawn of P. aeruginosa were performed. Note that D. discoideum actually feeds only on P. aeruginosa.

Table 2

Pseudomonas aeruginosa cells were harvested during exponential phase, resuspended into SM broth at a final O.D 600 of 1, diluted into 5 ml of SM broth and spread on SM agar plates. A total of 9 suspensions of D. discoideum cells in 5μl droplets corresponding to three-fold serial dilutions of the initial culture containing 4.0x10 6 cells/ml were then spotted on lawns of each P. aeruginosa strain. The most concentrated droplets therefore contained about 20.000 cells, whereas the most diluted ones contained less than 10 cells. The number of D. discoideum cells in each dilution was controlled using an automated cell counter (Countess, Invitrogen). Serial dilutions were performed in infection medium containing HL5 and 1x Soerensen buffer mix in 1:1 ratio [50 x Soerensen buffer, K 2 HPO 4 (99.86 g), Na 2 HPO 4 (17.8 g), add H 2 O 1000 ml]. The plates were incubated for 6 days at 19 • C to visualize the growth of D. discoideum.

Three patterns of Dictyostelium growth in the presence of P. aeruginosa could be easily recognized: inability to form phagocytic plaques on the bacterial lawn (pattern 0, Table 2), formation of phagocytic plaques without radial dissemination of cells (pattern 1), and formation of phagocytic plaques with radial dissemination of cells leading to the formation of an annular, transparent zone around the region were Dictyostelium cells were spotted (pattern 2). Only the most concentrated droplets containing 2 • 10 4 or 7 • 10 3 cells were able to form phagocytic plaques on the most virulent P. aeruginosa isolates (growth pattern 1, Table 2), whereas the presence of a transparent annular zone (growth pattern 2) was observed only in the less virulent P. aeruginosa isolates, therefore confirming the reliability and coherence of results obtained from both virulence assays. Dictyostelium cells were also spotted on a lawn of the non-pathogenic strain E. coli DH5α as a control. As expected, the amoeba was able to form phagocytic plaques on DH5α, even at the lowest cell concentrations and showed the formation of the annular zone (indicating radial dissemination) when 250 cells or more were spotted.

A mathematical model

The experimental results reported in Table 2 possibly reveal a complex dynamics underlying the interaction between amoeba and bacteria. In fact the occurrence of survival and extinction depending on virulence and on the initial concentration of amoeba suggests the existence of different states with roles changing as virulence changes.

The process we have to model, in order to describe the interaction amoebabacteria, does not fall within the traditional interactions between different species because, in our case, the amoeba is at the same time a predator (since it feeds on bacteria) and a victim (since bacteria belonging to a virulent strain are able to kill amoeba cells by their pathogenic action). Having in mind this main aspect, our model is based on the following description:

• in the absence of the amoeba, bacteria follow logistic growth, with intrinsic growth rate r and carrying capacity K;

• in the absence of bacteria, the amoeboid population undergoes an exponential decay with death rate m; in fact we suppose that bacteria are their unique source of food;

• amoeba cells feed on bacteria through a mass-action like mechanism with attack rate a;

• amoeba growth occurs at a rate proportional to the uptake of bacteria, with proportionality constant d;

• bacteria attack and kill amoeba cells with a functional response of Holling type with handling time T and attack rate b.

All the previous statements have a precise interpretation from the modeling point of view. Denoting respectively by u(t) and v(t) the number of bacteria and the number of amoeba cells, our model can be written as follows:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ u (t) = r 1 - u(t) K u(t) -au(t)v(t), v (t) = -mv(t) + du(t)v(t) - bu(t)v(t) 1 + bT v(t) .
(1)

By introducing the following non-dimensional variables

t = rt , ũ = u K , ṽ = a r v
and dropping all superscripts for simplicity, system (1) becomes

⎧ ⎪ ⎨ ⎪ ⎩ u (t) = (1 -u(t))u(t) -u(t)v(t), u (0) = u 0 , v (t) = -μv(t) + δu(t)v(t) -γ u(t)v(t) 1 + τv(t) , v(0) = v 0 , (2)
where

μ = m r , δ = dK r , γ = bK r , τ = bT r a
and u 0 and v 0 are the initial conditions.

This rescaling has reduced the number of parameters to four and we shall discuss the model dynamics with respect to γ representing the bacterial virulence. All other parameters will be considered as fixed. In particular, since our experimental framework suggests that amoeba life-time is very large if compared with the time scale of the interactions, we will assume μ 1.

(3)

Moreover since we have classified virulence on the basis of killing percentage in the first four hours (see Table 1), while the amoeba growth process is measured for 6 days, we will assume

τ 1 ( 4 )
to mean that the pathogenic process is rather fast.

A first analysis to determine the dynamics of the model concerns existence of steady states.

Steady states

Here we look for steady states of system (2), namely we look for nonnegative solutions of the system

⎧ ⎪ ⎨ ⎪ ⎩ (1 -u -v)u = 0, -μ + δu -γ u 1 + τv v = 0.
(5)

It is immediately seen that there are always the trivial steady state O ≡ (0, 0) and the bacteria-only state B ≡ (1, 0). However, the existence of coexistence steady states, for which amoeba and bacteria densities are both strictly positive, is not so straightforward and we need to analyze the system

⎧ ⎨ ⎩ 1 -u -v = 0, -μ + δu -γ u 1 + τv = 0, (6)
or, equivalently, the intersection of the curves

v = φ 1 (u) = 1 -u, v = φ 2 (u) = (γ -δ)u + μ τ (δu -μ) = 1 τ γu δu -μ -1 . (7)
Actually, we look for meaningful intersections of φ 1 (u) and φ 2 (u), i.e. such that

0 < u * < 1, 0 < v * < 1. (8)
Thus, noticing that

φ 2 (u) < 0 for 0< u < μ δ , φ 2 (u) > γ -δ τδ for u > μ δ ,
we conclude that

No coexistence equilibrium exists if at least one of the following conditions is satisfied

μ > δ (9) γ > δ(1 + τ ) (10)
However we note that we can exclude case (9) because condition (3) allows to restrict to the case

μ δτ 1 + τ < δ (11)
and, in view of (10), we will investigate existence of equilibria against different values of the parameter γ ∈ [0, δ(1 + τ )]. To this aim we transform system (6), by substitution, into

⎧ ⎨ ⎩ i) 1 -u -v = 0, ii) τδu 2 + [γ -δ(1 + τ ) -μτ ] u + μ(1 + τ ) = 0, ( 12
)
and solve equation (12, ii) to get

u ± (γ) = [δ(1 + τ ) -γ + μτ ] ± Δ(γ) 2δτ , (13)
where

Δ(γ) = [δ(1 + τ ) -γ + μτ ] 2 -4τδμ(1 + τ ).
We have indicated dependence on γ and we actually note that the discriminant Δ(γ) is a decreasing function of γ, on the interval [0, δ(1+τ )], vanishing at the point

γ + = δ(1 + τ ) - √ μτ 2 < δ(1 + τ ). ( 14
)
Figure 1 As a consequence, ( 13) provides two branches u ± (γ) of real positive solutions of (12, ii) for γ ∈ [0, γ + ], glued at γ + where

u -(γ + ) = u + (γ + ) = μ(1 + τ ) δτ < 1. ( 15
)
Moreover, u + (γ) is decreasing with γ and

u + (0) = 1 + 1 τ > 1, u + (γ * ) = 1,
where

γ * = δ -μ < γ + , (16)
so that u + (γ) satisfies ( 8) only for γ in the interval [γ * , γ + ]. Concerning the other branch u -(γ), we have that, since

u -(γ)u + (γ) = μ(1 + τ ), then u -(γ) is increasing and (see (15)) 0 < u -(γ) < u -(γ + ) < 1.
Thus u -(γ) satisfies ( 8) for γ ∈ [0, γ + ]. Finally, taking v ± (γ) = 1-u ∓ (γ), we have solutions of system ( 12) for γ ∈ [0, γ + ]. We note (see ( 12)) that v ± (γ)

are the two solutions of

δτ v * 2 + [δ -γ -τγ * ] v * + γ -γ * = 0 (17)
and that for γ < γ * only v + (γ) is positive.

In Figure 1 and Figure 2 we respectively show the bifurcation graphs of u ± (γ) and v ± (γ), with γ as a bifurcation parameter.

Figure 2

Next we discuss stability of the equilibria.

Stability

First of all we consider the trivial steady state O ≡ (0, 0) and the bacteriaonly state B ≡ (1, 0). The respective Jacobian matrices for these states are

J(O) = ⎛ ⎝ 1 0 0 -μ ⎞ ⎠ , J(B) = ⎛ ⎝ -1 -1 0 γ * -γ ⎞ ⎠ .

Consequently we have that

The trivial state O is always unstable, while B is stable if γ > γ * and unstable otherwise.

For analyzing the coexistence steady states we have to consider the Jacobian matrix

J(E * ) = ⎛ ⎜ ⎝ -u * -u * δ - γ 1 + τv * v * γτ u * v * (1 + τv * ) 2 ⎞ ⎟ ⎠ , (18)
and analyze the sign of trace and determinant of J(E * ) versus γ, along the different branches shown in Figure 2. Actually, we have

trJ(E * ) = u * (1 + τv * ) 2 τ 2 v * 2 + τ (2 -γ)v * + 1 det J(E * ) = u * v * (1 + T v * ) 2 δτ 2 v * 2 + 2δτ v * + δ -γ -γτ and, considering first the upper branch v + (γ) = 1 -u -(γ), for γ ∈ [0, γ + ],
we are led to analyze the sign of

T + (γ) = τ 2 v 2 + (γ) + τ (2 -γ)v + (γ) + 1 = = τ 2 (v + (γ) -t + (γ)) (v + (γ) -t -(γ)) , (19)
and

D + (γ) = δτ 2 v 2 + (γ) + 2δτ v + (γ) + δ -γ -γτ = = δτ 2 (v + (γ) -d + (γ)) (v + (γ) -d -(γ)) , (20)
where

t ± (γ) = γ -2 ± γ(γ -4) 2τ , d ± (γ) = -δ ± δγ(1 + τ ) δτ ,
looking for positive T + (γ) and D + (γ) in order to have stability.

First we prove that

D + (γ) is positive for γ ∈ [0, γ + ). (21)
In fact, first of all, for γ ≤ δ 1 + τ , d ± (γ) are both non-positive so that

D + (γ) = δτ 2 (v + (γ) + |d + (γ)|) (v + (γ) + |d -(γ)|) > 0.
To check the case γ > δ 1 + τ , we note that, since v -(γ) = 1u + (γ) and

v + (γ) = 1 -u -(γ)
are the roots of equation ( 17), then the quadratic polynomial

P (z) = δτ z 2 + [δ -γ -τγ * ] z + γ -γ *
can be written as

P (z) = δτ (z -v -(γ)) (z -v + (γ)) ,
so that, for z > 0,

P (z) < 0 if and only if v -(γ) < z < v + (γ).
Then, after some algebra and using ( 14), we get

τδP (d + (γ)) = 2δτ (1 + τ ) -(δ + γ + τγ * ) δγ(1 + τ ) = 2 δγ(1 + τ ) -γ -δ(1 + τ ) + τμ δγ(1 + τ ) = τμ - δ(1 + τ ) - √ γ 2 δγ(1 + τ ) < τμ - δ(1 + τ ) -γ + 2 δγ(1 + τ ) = 0, so that v -(γ) < d + (γ) < v + (γ)
and by ( 20), since d -(γ) < 0, we have (21).

Now we consider T + (γ) and prove

T + (γ) is positive for γ ∈ [0, γ + ). (22)
First we note that, for γ < 4, t ± (γ) are complex conjugate so that

4T + (γ) = (2τv + (γ) + 2 -γ) 2 + (4 -γ)γ > 0.
Thus, ( 22) is true if γ + ≤ 4, while if γ + > 4 we are left with the case 4 ≤ γ < γ + . To analyze this case, using

τ 2 t 2 ± (γ) + τ (2 -γ) t ± (γ) + 1 = 0,
we compute

P (t ± (γ)) = [(δ -1)γ -δ -τγ * ] t ± (γ) + γ -γ * - δ τ . (23)
Then, since τ < 1 (see ( 4)) and noticing that γ + ≥ 4 implies δ > 2, we have

P (t ± (4)) = 2μ + 2 τ (δ -2)(1 -τ ) > 0. ( 24
)
Now, it is easy to check that, for γ > 4, the functions t + (γ) and γt -(γ) are increasing with γ, while the function t -(γ) is decreasing. Moreover, since δ > 2, γ > 4 and τ < 1, we have

(δ -1)γ -δ -τγ * > (2δ -4) + δ(1 -τ ) + τμ > 0,
that, used in ( 23) implies that P (t ± (γ)) is increasing and, due to (24), it is positive for 4 ≤ γ < γ + . As a consequence, either

t -(γ) < t + (γ) < v -(γ) < v + (γ), or v + (γ) < t -(γ) < t + (γ)
In both cases T + (γ) > 0 and ( 22 As before we analyze the sign of

D -(γ) = δτ 2 v 2 -(γ) + 2δτ v -(γ) + δ -γ -γτ = = δτ 2 (v -(γ) -d + (γ)) (v -(γ) -d -(γ)) , (25)
with

d ± (γ) = -δ ± δγ(1 + τ ) δτ .

Figure 3

To discuss this point we represent the data of Tables 1 and2 in a more significant way introducing the variable

m = -ln (1 -k) ( 27
)
where k is the killing percentage of a strain as determined in Table 1. m is actually proportional to the mortality induced by the pathogenic strain and is then proportional to the parameter γ introduced in the model. In Figure 3 we represent the outcome of the culture at different values of m and of the number v 0 of amoeba cells. Namely, we label each point (m, v 0 ), available from Table 2, by an index 0, 1, 2, with the following meaning 0 = amoeba extinction , 1 = amoeba survival at a low density , 2 = amoeba survival at a high density .

( • for values greater than m = 4 amoeba does not survive whatever be the initial number of cells; this situation corresponds to values of γ greater than γ + , for which only the state B ≡ (1, 0) exists and is attractive.

The correspondence described in the present work is somewhat qualitative and could be better defined if the initial and final number of amoeba cells could be counted with more accuracy. This study still provides a first insight into the mathematical laws underlying the behavior of this complex biological system made of prey-predator and host-pathogen relationships.

The model and its bi-stability confirmed the importance of assessing the correct experimental parameters (mainly the initial number or concentration of amoeba cells) for an adequate application of the virulence assay for bacterial species or strains of unknown aggressiveness. Future work will focus on the

Figures Captions

• Caption for Figure 1 Bifurcation graph concerning existence of steady states. The component u(γ) corresponding to bacteria is represented as a function of γ.

Both branches u + (γ) and u -(γ), given in (13) are shown. γ * and γ + are respectively given in ( 16) and ( 14).

• γ * and γ + are respectively given in ( 16) and ( 14).

• Caption for Figure 3 Results from co-cultures of Dictyostelium amoebae with different strains of Pseudomonas aeruginosa. Different outcomes are represented according to (28) for different points (m, v 0 ) where m is defined in ( 27) and v 0 is the initial number of amoeba cells. After m = 5 the scale unit has been reduced.

  were performed with P. aeruginosa cells harvested during exponential phase, re-suspended into PG broth at a final O.D 600 of 1 and spotted on PG-agar plates in a final volume of 200μl on Peptone Glucose agar plates containing sorbitol at 0.15 M. The plates were initially incubated at room temperature for 18 hours. The adult nematodes (on average 60 -70 individuals) were grown on E. coli OP50 on NGM plates, picked and placed on the P. aeruginosa spots. As a control, one plate was seeded with the non-pathogenic E.coli strain DH5α. The plates were subsequently incubated at 22 • C for 4 hours. The viability/mortality of C. elegans was observed under an optical microscope at 10X magnitude. Non-motile individuals were considered as dead, and the killing percentage was calculated as the ratio between dead and total C. elegans individuals. Clinical isolates of P. aeruginosa showed a wide range of virulence against C. elegans.

  ) is true. Now we consider the lower branch of Figure 2, v -(γ) = 1-u + (γ), existing only for γ ∈ [γ * , γ + ].

)

  Figure 3 we see that

  Caption for Figure 2 Bifurcation graph concerning existence of steady states. The component v(γ) = 1-u(γ) corresponding to the amoeba is represented versus γ. Both branches v + (γ) = 1u -(γ) and v -(γ) = 1u + (γ) are shown.

Table 1 :

 1 C. elegans killing assays: different strains of P. aeruginosa tested agains C. elegans. The killing percentage is measured after 4 hours.

	strain	% killing
	Hpu45	100,00
	N-10	98,99
	N-54	98,83
	N5	96,77
	N1	95,30
	Hpu49	95,30
	N-89	94,15
	N101	94,00
	N47	93,30
	Hpu26	92,50
	N-3	90,65
	N60	89,50
	N-36	89,20
	N-81	85,76
	Hpu75	84,84
	Hpu105	79,94
	Hpu47	78,00
	N-85	76,00
	Hpu92	70,6
	N106	68,61
	Hpu28A	61,30
	N-30	52,64
	N-11	45,47
	N29	29,00
	Hpu28	27,60
	Hpu56	26,96
	N-116	6,66
	DH5-α	0,00

1 5. Tables

✩ This work was supported by MIUR (Ministero Istruzione Università Ricerca) within the PRIN 2007 project "Mathematical Population Theory: methods, models,

The previous results can be summarized in the following statement on the coexistence state

The coexistence state

is stable, while the state

is unstable.

Comparing experiments and model behavior

The previous analysis of the model reveals that the so called bi-stable behavior, for a certain range of the parameter representing virulence, is a characteristic feature of the amoeba-bacteria interaction. A qualitative comparison of the model with the experimental results shows that the behavior is actually confirmed by data.

assessment of growth properties of bacteria and amoeba separately, including spatial patterns.

Acknowledgements

The authors wish to thank Andrea Pugliese for his many suggestions and comments during the preparation of the manuscript.