Viola Pavlová
email: viola.pavlo@gmail.com

Ludék Berec
email: berec@entu.cas.cz

David S Boukal
email: boukal@entu.cas.cz

Luděk Berec

Luděk Berec

Caught between two allee effects: Trade-off between reproduction and predation risk

Keywords: anti-predator behaviour, dilution effect, extinction, mate search, population resilience

Reproductive activities are often associated with conspicuous morphology or behaviour that could be exploited by predators. Individuals can therefore face a trade-off between reproduction and predation risk. Here we use simple models to explore population-dynamical consequences of such a trade-off for populations subject to a mate-finding Allee effect and an Allee effect due to predation. We present our results in the light of populations that belong to endangered species or pests and study their viability and resilience. We distinguish several qualitative scenarios characterized by the shape and strength of the trade-off and, in particular, identify conditions for which the populations survive or go extinct. Reproduction can be so costly that the population always goes extinct. In other cases, the population goes extinct only over a certain range of low, intermediate or high levels of reproductive activities. Moreover, we show that predator removal (e.g. in an attempt to save an endangered prey species) has the least effect on populations with low cost of reproduction in terms of predation and, conversely, predator addition (e.g. to eradicate a pest) is most effective for populations with high predation cost of reproduction. Our results indicate that a detailed knowledge of the trade-off can be crucial in applications: for some trade-off shapes, only intermediate levels of reproductive activities might guarantee population survival, while they can lead to extinction for others. We therefore suggest that the fate of populations subject to the two antagonistic Allee effects should be evaluated on a case-by-case basis. Although the literature offers no quantitative data on possible trade-off shapes in any taxa, indirect evidence suggests that the trade-off and both Allee effects can occur simultaneously, e.g. in the golden egg bug Phyllomorpha laciniata.

Introduction

Various trade-offs affect individual life history decisions and in turn influence dynamics of populations. One such trade-off links reproduction and predation risk.

Many species have evolved anti-predator behaviour during which they stop mating and perform an escape manoeuvre (e.g. [START_REF] Svensson | The odour makes the difference: male moths attracted by sex pheromones ignore the threat by predatory bats[END_REF][START_REF] Svensson | Listening in pheromone plumes: disruption of olfactory-guided mate attraction in a moth by a bat-like ultrasound[END_REF]. Higher levels and/or prolonged periods of spatial movement [START_REF] Anholt | Interaction between food availability and mortality mediated by adaptive behaviour[END_REF][START_REF] Kotiaho | Male mating success and risk of predation in a wolf spider: balance between sexual and natural selection?[END_REF] or sexual signalling [START_REF] Zuk | Exploitation of sexual signals by predators and parasitoids[END_REF] during mate search help find a mate or choose a better one but can also attract predators. For example, neotropical bats can exploit frog calls [START_REF] Ryan | Bat predation and sexual advertisement in a neotropical anuran[END_REF][START_REF] Tuttle | Bat predation and the evolution of frog vocalizations in the neotropics[END_REF], parasitic flies can target calling male field crickets [START_REF] Lewkiewicz | Latency to resume calling after disturbance in the field cricket, Teleogryllus oceanicus, corresponds to population-level differences in parasitism risk[END_REF] and skuas can respond to male petrel calls [START_REF] Mougeot | Predation as a cost of sexual communication in nocturnal seabirds: an experimental approach using acoustic signals[END_REF]. Predation risk may also be high during copulation [START_REF] Trochine | When prey mating increases predation risk: the relationship between the flatworm Mesostoma ehrenbergii and the copepod Boeckella gracilis[END_REF], pregnancy, spawning and breeding period when the individuals are often less motile or easier to detect. For example, gravid females of the copepod Eudiaptomus gracilis become an easier prey as the brood looks like a conspicuous dark spot inside the body [START_REF] Svensson | Fish predation on Eudiaptomus gracilis in relation to clutch size, body size, and sex: a field experiment[END_REF] and gravid females of other copepods can be eaten by fish because of reduced manoeuvrability [START_REF] Winfield | The cost of copepod reproduction: increased susceptibility to fish predation[END_REF].

Both reproductive success and probability of avoiding predation, involved in the reproduction-predation risk trade-off, may be positively related to population size or density and therefore subject to a component Allee effect [START_REF] Courchamp | Inverse density dependence and the Allee effect[END_REF][START_REF] Courchamp | Allee effects in ecology and conservation[END_REF]Stephens and Sutherland, 1999). Allee effects are mostly studied as single mechanisms of positive density dependence, but multiple Allee effects can be common and the ways in which two or more component Allee effects interact can be complex [START_REF] Berec | Multiple Allee effects and population management[END_REF].

Various mechanisms can lead to component Allee effects in reproduction (Table 1).

In particular, individuals in small or sparse populations can face a mate-finding Allee effect and have reduced reproductive success due to their inability to find a (suitable) mate during their receptive period. The mechanisms include male-female encounters in motile animals, encounters between pollen grains and egg cells in plants, and encounters between male and female gametes in many sessile or semi-sessile marine A c c e p t e d m a n u s c r i p t 4 invertebrates [START_REF] Gascoigne | Dangerously few liaisons: a review of mate-finding Allee effects[END_REF]. Evidence for Allee effects affecting reproduction in general and for mate-finding Allee effects in particular is quite widespread; observations of the latter span a wide spectrum of taxa from sheep ticks [START_REF] Rohlf | The effect of clumped distributions in sparse populations[END_REF] to polar bears (Molnár et al., 2008) among terrestrial animals and from copepods [START_REF] Kiørboe | Sex, sex-ratios, and the dynamics of pelagic copepod populations[END_REF] to queen conch [START_REF] Stoner | Evidence for Allee effects in an over-harvested marine gastropod: density-dependent mating and egg production[END_REF] among marine species. Other mechanisms by which larger population size or density may enhance reproduction success include sperm limitation, physiological stimulation of reproduction and female choice [START_REF] Gascoigne | Dangerously few liaisons: a review of mate-finding Allee effects[END_REF].

Likewise, individuals in small or sparse populations may be more susceptible to predation due to reduced ability to avoid or escape predators thanks to passive or active defence. Prey can actively defend themselves through predator mobbing [START_REF] Krebs | An introduction to behavioural ecology[END_REF], more vigilant scanning for predators [START_REF] Mooring | Vigilance, predation risk and the Allee effect in desert bighorn sheep[END_REF] or confusion effect when schooling [START_REF] Nottestad | Herring schooling manoeuvres in response to killer whale attacks[END_REF]. Passive defence (the dilution effect) is characterized by per-capita predation probability declining with increasing prey density and corresponds to a type II functional response [START_REF] Gascoigne | Allee effects driven by predation[END_REF]. For example, woodland caribou face higher predation risk in smaller groups as predators usually kill one animal per group per attack irrespective of the group size [START_REF] Wittmer | The role of predation in the decline and extirpation of woodland caribou[END_REF].

In this article, we explore the dynamics of a population that faces a trade-off between two Allee effects: one due to reproduction and the other due to predation. Our study is motivated by situations in which the population belongs to an endangered species, and we identify properties of its predators and trade-offs that might put the population under extinction risk. Alternatively, our model describes a pest species that needs to be controlled or eradicated. To quantify the extinction risk and eradication potential we use population resilience, defined as the maximum disturbance the population in a stable state may sustain to avoid extinction [START_REF] Beisner | Alternative stable states in ecology[END_REF], and relative population resilience, by which we mean a relative change in the resilience of the population after a predator removal/addition. When measuring relative resilience, we also distinguish prey with flexible and inflexible reproductive behaviour. We assume that flexible reproductive behaviour (determining strength of the mate-finding Allee effect) can be instantly changed after predator removal or addition, whereas inflexible behaviour cannot respond to predator presence or absence. They represent extreme

A c c e p t e d m a n u s c r i p t

5 but useful approximations of real reproductive behaviour, which can involve both morphological and behavioural traits. Most morphological traits such as bright colouration of males in many birds and many types of behaviour such as lekking displays of various birds and flies [START_REF] Andersson | Sexual selection[END_REF] are 'hardwired' characteristics that individuals cannot change and therefore correspond to the inflexible behaviour.

On the other hand, quite a number of behavioural traits such as mating calls in some orthopteroid insects are plastic. They can be adjusted to the perceived predation risk [START_REF] Zuk | Exploitation of sexual signals by predators and parasitoids[END_REF] and we represent them by the flexible behaviour in our model.

We ask the following questions: When do the two component Allee effects lead to extinction of the population and how do they affect its extinction risk or eradication potential? How are the results affected by the shape of the reproduction-predation risk trade-off? And finally, are the results for prey with flexible and inflexible reproductive behaviour different?

Methods

We use several simplifying assumptions on the life histories and population densities of males and females in this paper: identical mortalities, balanced sex ratio at birth and equal initial densities. This allows us to follow changes of the total population density N in time without discerning between males and females [START_REF] Boukal | Modelling mate-finding Allee effects and populations dynamics, with applications in pest control[END_REF]:

P N F K N dN N M bN dt dN ) ( 1 ) ( - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + - = (1)
That is, we assume that negative density dependence affects only the prey mortality rate but not its birth rate and that predators do not interfere with each other and do not 2.

Allee effects in reproduction and due to predation

Mate-finding Allee effects are commonly quantified through positive density dependence in the female mating rate [START_REF] Boukal | Single-species models of the Allee effect: extinction boundaries, sex ratios and mate encounters[END_REF][START_REF] Courchamp | Allee effects in ecology and conservation[END_REF][START_REF] Dennis | Allee effects: population growth, critical density and the chance of extinction[END_REF]. The simplest and at the same time most widely used models of positively density-dependent mating rates include the exponential function M(N) = 1exp(-N/θ) and the hyperbolic function M(N) = N/(N+θ), in which θ quantifies the strength of mate-finding Allee effect [START_REF] Dennis | Allee effects: population growth, critical density and the chance of extinction[END_REF]. We use the hyperbolic function in model ( 1) as it allows for some analytical calculations. Any of these functions can likewise be used as a (phenomenological) model of other Allee effects in reproduction beyond mate finding [START_REF] Courchamp | Allee effects in ecology and conservation[END_REF].

Predation-driven Allee effects arise for type II functional responses. We use two alternative descriptions, F(N) = αN/(1+N/β) and F(N) = λN/(1+ λhN). The latter form is the standard formulation due to [START_REF] Holling | Some characteristics of simple types of predation and parasitism[END_REF], in which λ scales the predator-prey encounter rate and h is the handling time of one prey individual, and represents an Allee effect due to predator dilution [START_REF] Berec | Multiple Allee effects and population management[END_REF][START_REF] Courchamp | Allee effects in ecology and conservation[END_REF][START_REF] Gascoigne | Allee effects driven by predation[END_REF]. The former expression can be interpreted as a predation-driven Allee effect invoked by a behavioural response of the population (e.g. herding or mobbing) to a predator with an otherwise linear functional response:

the behavioural response is more efficient at higher population densities and reduces the (linear) attack rate α by a factor of 1/(1 + N/β). Note that β scales the maximum per-capita attack rate (= αβ) but also defines the population density at which F(N)

declines to 50% of its maximum value, and we call β prey vulnerability in this article.

Inserting functional forms for M(N) and F(N) to model (1), our alternative models of population dynamics become:

β α θ / 1 1 N NP K N dN N N bN dt dN + - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + - + = (2a)
+ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + - + = 1 1 (2b)
We rescale both models to reduce the number of parameters from seven to five. In model (2a), we rescale the state variable as x = N/K and time as τ = td and define new parameters Θ = θ/K (relative strength of mate-finding Allee effect), B = β/K (relative prey vulnerability), g = b/d and δ = P/d, to obtain

B x x x x x x gx d dx / 1 ) 1 ( + - + - Θ + = δ α τ (3a)
Rescaling model (2b) analogously, we get

Hx x x x x x gx d dx λ δ λ τ + - + - Θ + = 1 ) 1 ( (3b)
where H = hK (relative handling time) and the other parameters are as above. The parameters α, B and λ scale the strength of Allee effect due to predation: predation mortality in models (3a) and (3b) increases with increasing α, B and λ. Only the parameter Θ scales the strength of mate-finding Allee effect in these models, and reproductive success declines with increasing Θ. We skip the adjective 'relative' when further referring to parameters Θ, B and H.

Reproduction-predation risk trade-off

Trade-off shapes can strongly influence evolutionary dynamics and endpoints [START_REF] Bell | The costs of reproduction and their consequences[END_REF][START_REF] Rueffler | Adaptive walks on changing landscapes: Levins' approach extended[END_REF][START_REF] Hoyle | The influence of trade-off shape on evolutionary behaviour in classical ecological scenarios[END_REF]. We thus consider a range of trade-off shapes that cover varying costs of reproduction in terms of predation risk, and investigate two trade-offs for model (3a),

A c c e p t e d m

a n u s c r i p t

8 z z / 1 max max 1 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Θ Θ - = α α (4a) and z z B B / 1 max max 1 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Θ Θ - = , (4b)
and one trade-off for model (3b),

z z / 1 max max 1 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Θ Θ - = λ λ . ( 4c

Model analysis

Models (3a) and (3b) can have, in addition to the extinction equilibrium x 0 = 0, either no positive equilibrium or a unique positive equilibrium that is stable or two positive equilibria, of which the higher (labelled S) is stable and the lower (U) unstable (see below for exact results). We focus here only on the difference R = S-U and call it population resilience. It can be interpreted as the maximum disturbance after which the population can still return back to the stable equilibrium without going extinct. If there is only a unique positive equilibrium that is stable, we set U = 0.

We start with standard viability analysis of models (3a) and (3b) and search for the stable and unstable equilibria and conditions for their existence. In the simplest cases when Θ and/or α or λ are zero or B infinitely large (that is, one or both component Allee effects are absent), the equilibria can be calculated analytically. We use MATLAB (The MathWorks, Inc.) to find the equilibria numerically when both Allee effects are present.

To investigate the role of reproduction-predation risk trade-off, we examine the impact of changes in the strength of mate-finding Allee effect Θ on the population resilience R along the trade-off in the AR, PV and ER models. In addition, we explore relative population resilience along the trade-off for prey with flexible and inflexible reproductive behaviour. In line with the general considerations above, by 'flexible'

we mean that prey reproductive behaviour can be changed instantly such that there is no mate-finding Allee effect (Θ = 0) in the absence of predators and there is some (Θ > 0) in their presence, while prey with inflexible reproductive behaviour are assumed to keep the same Θ > 0 regardless of predators. For prey with flexible behaviour, we therefore measure the relative population resilience as ρ f = (S-U)/(S 0 -U 0 ), where S 0 and U 0 correspond to no predation (P = 0) and no mate-finding Allee effect (Θ = 0), i.e. U 0 = 0; in fact ρ f is just a multiple of the population resilience R. For prey with inflexible behaviour, we measure the relative population resilience as ρ i = (S-U)/(S 1 -

A c c e p t e d m a n u s c r i p t

10 U 1 ), where S 1 and U 1 correspond to no predation and the mate-finding Allee effect remaining at the unchanged strength Θ > 0; ρ i always equals 1 if there is no predation.

Results

Equilibria of models (3a) and (3b) and their stability

The extinction equilibrium x 0 = 0 of model ( 3a) is locally stable for any mate-finding Allee effect (Θ > 0). It is globally stable, i.e. the population goes extinct regardless of its density, and no other equilibria exist if the mate-finding Allee effect and/or the predation-driven Allee effect are sufficiently strong (i.e. if Θ and/or α are sufficiently large). For no mate-finding Allee effect (Θ = 0), x 0 is locally or globally stable if α > (g-1)/δ and unstable if α < (g-1)/δ. In the latter case, the population always reaches a unique positive equilibrium (environmental carrying capacity) regardless of its initial density. A complete analysis of the case with no matefinding Allee effect is given in Appendix (Table A1).

If the extinction equilibrium x 0 = 0 is only locally stable, two positive equilibria exist of which the lower is unstable (Allee threshold) and the upper is locally stable (carrying capacity). All else being equal, the two positive equilibria approach one another with increasing Θ and/or α until both merge and the extinction equilibrium becomes globally stable (Fig. 2A). The population resilience R is maximized for no mate-finding Allee effect and no predation (Θ = 0 and α = 0; R = g-1) and declines as the strength of mate-finding Allee effect Θ and/or the attack rate α increase (Fig. 3A;

the thin lines depict isolines along which the resilience is constant). The lower limit of the isolines is the viability limit at which the positive equilibria merge and R = 0. For Θ and/or α to the right and above the viability limit, the population cannot persist as deaths exceed births for any population density. We denote the mate-finding Allee effect strength and the attack rate at which the viability limit intersects the respective axes as Θ lim and α lim (Θ lim = 0.172 and α lim = 2.55 in Fig. 3A). The relative population resilience ρ i of prey with inflexible reproductive behaviour is maximized in the absence of predation (α = 0) and declines with increasing α for any 0 < Θ < Θ lim (Fig. 3B).

A c c e p t e d m a n u s c r i p t

11 For high attack rates, α > (g-1)/δ, the dependence of equilibria of model (3a) on the mate-finding Allee effect strength Θ and the prey vulnerability B is qualitatively the same as in Fig. 2A and the corresponding (relative) population resilience plots are analogous to Fig. 3A-B. For low attack rates, α < (g-1)/δ, the results differ qualitatively and high values of B do not affect the equilibria (Fig. 2B). This is because for highly vulnerable prey (1/B close to 0) the predation term in model (3a) reduces to αδx < (g-1)x, and hence for Θ = 0 the population always attains a unique positive equilibrium (carrying capacity) that is globally stable. In biological terms, the prey cannot go extinct if the mate-finding Allee effect is absent and the attack rate is low enough. Consequently, the viability limit and isolines of the (relative) population resilience in the Θ-B parameter space are much more convex (Fig. 3C-D).

Similar conclusions hold for the existence and stability of the equilibria of model (3b).

The extinction equilibrium x 0 = 0 is locally or globally stable for any mate-finding Allee effect (Θ > 0). For no mate-finding Allee effect (Θ = 0), x 0 = 0 is locally or globally stable for high prey encounter rates λ > (g-1)/δ and unstable for λ < (g-1)/δ. A complete analysis for Θ = 0 (see Table A2 in Appendix) implies that if the handling time is sufficiently short, H < 4δ /(g-1) 2 , the positive equilibria of model (3b), population resilience and relative population resilience behave similarly as shown in Figs 2A, 3A and 3B, respectively. On the other hand, long handling times, H > 4δ /(g-1) 2 , yield results that are qualitatively analogous to Figs 2B, 3C and 3D and the encounter rate λ affects the location of equilibria only when relatively small. This is because the predation-driven Allee effect is not strong enough to drive the population to inevitable extinction even for virtually infinite encounter rates λ, for which the predation term in model (3b) reduces to a constant δ/H.

The isoline plots of (relative) population resilience in Fig. 3 also reveal that the prey population can or cannot always be eliminated by strong predation. Populations that can always be eliminated by strong predation are characterized by plots in which the viability limit intersects the y-axis at some α lim < ∞ (Fig. 3A-B). These populations will be wiped out by predation with α > α lim irrespective of the strength of matefinding Allee effect Θ; we call them populations with predation limit. On the other A c c e p t e d m a n u s c r i p t 12 hand, the viability limit of some populations may not intersect the y-axis (Fig. 3C-D).

These populations cannot be wiped out by predation of any strength when the matefinding Allee effect strength Θ is sufficiently low and we call them populations without predation limit. Populations without predation limit appear in the PV model with low attack rates α and in the ER model with long handling times H (Table 3).

Population resilience under reproduction-predation risk trade-off

How do population resilience and relative population resilience change along the reproduction-predation risk trade-off? And how is the population affected by the shape of the trade-off? We first explain several examples, based on the AR model, in detail before giving the complete results. We have already shown how the population resilience changes with the strength of mate-finding Allee effect and the predation pressure (Fig. 3A). We now overlay different trade-off shapes over the isoline landscape. The curves A, B and C in Fig. 4A are three trade-off curves (4a) that exemplify the main differences in the qualitative outcomes reported in this article.

First, there is a difference between cheap and costly reproduction. For sufficiently 'cheap' reproductive activities characterized by a convex trade-off curve (z < 1, curve A in Fig. 4A), the resilience R can be maximized at intermediate values of the matefinding Allee effect strength (curve A in Fig. 4B). Populations with more costly reproductive activities characterized by linear and concave trade-offs (z ≥ 1, curve B in Fig. 4A) have the resilience always maximized either in the absence of predation (Θ = Θ max and α = 0) or in the absence of the mate-finding Allee effect (Θ = 0 and α = α max ; as for curve B in Fig. 4B). The exact value of z above which the maximum resilience can no longer occur at an intermediate Θ depends on the curvature of the trade-off relative to the curvature of the isolines; the value will be z < 1 because the isolines are always convex.

Second, populations might be always safe or liable to extinction as the strength of mate-finding Allee effect and the predation pressure vary along the trade-off. This is determined by the location of the trade-off relative to the viability limit: populations characterized by trade-off curves A and B persist for any Θ between 0 and Θ max and 4B; the maximum strengths of the two Allee effects due to mate-finding and predation are larger for curve C than for curve B). In extreme cases, the trade-off can be so severe that the population always goes extinct (see below).

Finally, we note that predation ceases and only the mate-finding Allee effect operates for Θ > Θ max . Any two trade-offs thus give the same results for Θ above the larger of their Θ max values. In addition, if Θ max < Θ lim , the population resilience declines with increasing Θ for Θ > Θ max until it reaches zero at Θ = Θ lim (all curves in Fig. 4B).

Now we give the complete results. Results for populations with predation limit can be classified into five qualitatively different scenarios, characterized by the location of the point corresponding to the maximum double Allee effect without any trade-off, i.e. point [Θ max , α max ] or [Θ max , B max ] for model (3a) and point [Θ max , λ max ] for model (3b), relative to the viability limit. The five scenarios correspond to the five areas marked as A to E in Fig. 5A. Results for populations without predation limit can be described by only three scenarios (areas A to C in Fig. 5B), which are covered by the former five: parameter combinations from areas A to C in Fig. 5A and Fig. 5B give qualitatively analogous results. We thus present and discuss only the results for populations with predation limit.

In all five scenarios, the population resilience declines with increasing cost of reproduction (i.e. increasing z) for the mate-finding Allee effects satisfying Θ < Θ lim (Fig. 6). The main distinguishing feature of each scenario is whether and for what values of Θ the population cannot persist. The population will survive predation irrespective of the strength of mate-finding Allee effect and the shape of the trade-off only in scenario A (Fig. 6A). In biological terms, this scenario corresponds to prey that can survive predation even when both Allee effects are at maximum strength. In all other scenarios, the population may go extinct; the range of Θ that lead to extinction increases with z. In scenario B, the population will not survive predation for intermediate Θ and sufficiently costly reproduction (i.e. concave trade-offs with high values of z; Fig. 6B). This contrasts with the remaining three scenarios in which E, respectively (Fig. 6C-E). In other words, extinction occurs for strong mate-finding Allee effects (high Θ) if their maximum strength is very high (cases C and E) and for weak mate-finding Allee effects (low but positive Θ) if the maximum strength of Allee effect due to predation is very high (cases D and E). In the latter case, the tradeoff results in strong Allee effects due to predation. Finally, the population is never viable for sufficiently costly reproduction in scenario E (Fig. 6E).

The relative population resilience ρ f of prey with flexible reproductive behaviour is a multiple of the population resilience R and can thus be taken from Fig. 6. We observe the same five qualitatively different scenarios for the relative population resilience ρ i of prey with inflexible reproductive behaviour (Fig. 7), although the value of Θ at which the relative population resilience ρ i is maximized (or minimized) can be different from that for the (relative) population resilience R (ρ f ) (compare the corresponding panels in Figs. 6 and7).

Discussion

A trade-off between reproduction and predation risk occurs in many animals but its impact on population viability has been little studied. We have focused in this article on the viability of a population that faces, along with this trade-off, two different Allee effects due to reproduction and predation, and examined how various characteristics of the population and shapes of the trade-off affect its resilience to disturbances. We have also compared the results for prey that can or cannot adjust their reproductive activity to the predation level.

Summary of outcomes and their implications for population management

We have found that the fate of a population, given its demographic rates g and δ,

depends on three main factors: the cost of reproduction (that is, the trade-off shape characterized by the parameter z), the maximum strength of mate-finding Allee effect at which predation ceases and the maximum strength of Allee effect due to predation

A c c e p t e d m a n u s c r i p t

15 in the absence of the mate-finding Allee effect (that is, the trade-off location characterized by the parameters α max , B max or λ max , respectively). While the results depend strongly on the reproduction-predation risk trade-off, the flexibility or inflexibility of prey reproductive behaviour as well as the type of Allee effect due to predation (= models AR, PV and ER) do not produce qualitatively different outcomes.

We summarize all results jointly in what follows.

Overall, the (relative) population resilience declines with increasing cost of reproduction embodied in the parameter z that represents the strength of reproductionpredation risk trade-off. In some cases, reproduction can be so costly in terms of predation that the population always goes extinct. In other cases, the population goes extinct only over a certain range of mate-finding Allee effects, which may correspond to high, low or intermediate levels of reproductive activity. For example, intermediate levels of reproductive activity might be the only ones that guarantee population survival (e.g., sufficiently convex trade-offs in scenario E in Fig. 6E) or that, on the contrary, lead to population extinction (e.g., sufficiently concave trade-offs in scenario B in Fig. 6B).

Our results on the relative population resilience imply that the impact of predator manipulation (i.e., removal or addition) differs quantitatively but not qualitatively for prey with flexible and inflexible reproductive behaviour. That is, our main results are equally valid if predators exploit hardwired traits associated with reproduction (usually colouration or canalized behaviour) or if they exploit adjustable prey behaviour such as the intensity and duration of mating calls. All else being equal, predator manipulation has the largest effect on populations with high cost of reproduction in terms of predation. This can have either positive or negative consequences depending on the application: conservation biologists have obviously opposing goals to pest managers.

Predator removal is often considered among conservation management strategies because it can increase survivorship in endangered populations [START_REF] Sinclair | Predicting effects of predation on conservation of endangered prey[END_REF].

Our modelling results show that removal of predators can be particularly worthwhile for populations with predation limit that can be heavily predated and must suppress

A c c e p t e d m a n u s c r i p t

16 their reproductive behaviour considerably in order to avoid predation (i.e. have relatively high α max and Θ max ; area E in Fig. 5A). Such populations cannot survive without removal of the predators. On the contrary, pest management might attempt to remove a pest by introducing additional or entirely new predators [START_REF] Solomon | Biocontrol of pests of apples and pears in northern and central Europe -3. Predators[END_REF]. Here we demonstrate that, within the types of life histories and ecological interactions considered in our models, the most difficult pests fall into the category of populations that are without predation limit and face relatively few difficulties in finding mates when they manage to avoid predation (i.e. have relatively small Θ max ;

area A in Fig. 5B). Such pests are impossible to eradicate by adding predators.

Given the disparate results on population viability that depend on quantitative details of the reproduction-predation risk trade-off, it is clear that no general rule of thumb can be put forward and the fate of a population subject to antagonistic mate-finding and predation-driven Allee effects should be evaluated on a case-by-case basis.

Evidence for the main model ingredients

Although our results show that the shape and location of the reproduction-predation risk trade-off can be the key for survival or extinction of the population, we have found that quantitative data on it are virtually lacking in the literature. In fact, nearly all evidence supporting the existence of such a trade-off is based on 'binary' experiments where the prey reproductive activity is observed under the presence or absence of its predators (see the examples in Introduction). This indicates that the relationship is either difficult to measure or that it has so far largely escaped an interest of behavioural ecologists. Nevertheless, many species across various taxa (predominantly insects, but also crabs, fish, frogs or birds) exhibit a trade-off between reproduction and survival from predation, suggesting that the occurrence of the tradeoff might be common.

In addition to the presence of a reproduction-predation risk trade-off, the application of our model requires that the species be simultaneously subject to a mate-finding Allee effect and an Allee effect due to predation. [START_REF] Kramer | The evidence for Allee effects[END_REF] reviewed the empirical evidence for Allee effects and provided an extensive list of species likely subject to the Allee effect due to predation satiation and an Allee effect due to mate Our study might appear irrelevant in the light of this heterogeneous body of evidence.

However, as [START_REF] Kramer | The evidence for Allee effects[END_REF] suggest and we fully agree, expectations and studies of a certain type of the Allee effect can be biased in various species or groups of species (e.g. predation satiation is mostly proposed for herbivorous mammals and mate finding for invertebrates). Studies dealing with 'less obvious' Allee effects and multiple Allee effects among all species are lacking. Another reason might be that some combinations of taxa and Allee effects are more amenable to experimental manipulation or field observations than others. For example, insects are good model organisms for the mate-finding Allee effect, while the impact of predation and potential Allee effects driven by it are reasonably easily measured in small populations of large ungulates.

Simultaneous occurrence of both types of Allee effect has been put forth in [START_REF] Berec | Multiple Allee effects and population management[END_REF]. The candidate species include the monarch butterfly Danaus plexippus, in which individuals in small overwintering groups suffer mating depression in spring and less efficient predator dilution; the social spider Anelosimus eximius, where small colonies have a reduced number of eggs per sac and reduced offspring survival to maturation, owing to reduced protection against predators; and the white abalone Potential occurrence of all three mechanisms in a species can also be inferred indirectly from known details on its life history. For example, the golden egg bug

A c c e p t e d m a n u s c r i p t 18

Phyllomorpha laciniata exhibits a component Allee effect in reproduction: females produce twice as many eggs if there are enough conspecifics to lay eggs upon their backs [START_REF] Garcia-Gonzales | Oviposition site selection and oviposition simulation by conspecifics in the golden egg bug (Phyllomorpha laciniata): implications for female fitness[END_REF][START_REF] Katvala | The effect of conspecific density on female reproduction in an egg-carrying bug[END_REF]Kaitala, 2003, 2005) and the eggs survive better if carried on backs than if laid on leaves [START_REF] Garcia-Gonzales | Oviposition site selection and oviposition simulation by conspecifics in the golden egg bug (Phyllomorpha laciniata): implications for female fitness[END_REF]. Major predators of this species are ants, some of which are generalist predators [START_REF] Hölldobler | The ants[END_REF] and we might expect that they have a saturating functional response. Moreover, individuals that carry eggs are more visible and in a greater risk of predation [START_REF] Kaitala | Egg load and mating status of the golden egg bug affect predation risk[END_REF], copulation with an individual carrying eggs increases predation risk for both partners [START_REF] Kaitala | Egg carrying attracts enemies in a cryptic coreid bug (Phyllomorpha laciniata)[END_REF] and the frequency of copulations in the presence of a predator is reduced [START_REF] Miettinen | The effect of ant predation on survival and mating frequency of the golden egg bug in field experiments (Hymenoptera:Formicidae; Heteroptera:Coreidae)[END_REF]. Finally, copulation lasts for many hours and pairs are less mobile [START_REF] Miettinen | Do egg carrying and protracted copulation affect mobility in the golden egg bug?[END_REF] and very much susceptible to predation [START_REF] Kaitala | Egg load and mating status of the golden egg bug affect predation risk[END_REF]. We might therefore conclude that all three mechanisms probably occur in the golden egg bug; however, the available quantitative data are insufficient to apply our model to this species.

On the other hand, other adaptations to increased predation mortality might counteract the decrease in reproductive activities under predation risk. For example, a number of traits related to mating, such as the courtship rate, intensity of colour patterns or female preference for brightly coloured males, decline in response to increased predation risk in the guppy Poecilia reticulata. However, guppies from high-predation sites also invest relatively more energy to reproduction, have shorter interbrood intervals and produce more (and smaller) offspring compared to fish from lowpredation sites (Magnhagen, 1991, and references therein). Such additional responses to increased predation risk go beyond the reproduction-predation risk trade-off as understood in this article.

Concluding remarks

As our results suggest, searching for multiple Allee effects and their relationship in plant and animal populations could reveal new and interesting insights and is worth further research both by theoreticians and field workers. We foresee two distinct ways in which this theory can be extended. First, one could critically assess some of the simplifying assumptions that we made in the model. One restriction of our model is that all predation unrelated to reproductive activities of the population can be included in the mortality parameter d. We also assumed identical male and female life histories to simplify the analyses. Although both assumptions are often reasonable, there are many examples of sex-selective predation [START_REF] Magnhagen | Predation as a cost of reproduction[END_REF][START_REF] Zuk | Exploitation of sexual signals by predators and parasitoids[END_REF] and such predation can fundamentally change the predator-prey dynamics [START_REF] Boukal | Does sex-selective predation stabilize or destabilize predator-prey dynamics?[END_REF]. Second, one could examine evolutionarily stable reproductive strategies in the present setting to see how individuals benefit from putting energy to reproduction or to predator avoidance, a topic on which we currently work.

Last but not least, we have also highlighted how surprisingly little is known on the processes that, according to our results, drive population dynamics under the reproduction-predation risk trade-off. More data are needed to assess which shapes of the trade-off are most common in nature and how often species are actually caught between two Allee effects caused by reproduction and predation. Setting the right hand side of equation (3a) to zero always has the trivial solution x 0 = 0, corresponding to population extinction. The extinction equilibrium x 0 is locally stable for Θ > 0 as the derivative of the right hand side of equation (3a) with respect to x, evaluated at x 0 , is -1-αδ < 0 (all model parameters are positive).

To find positive model equilibria for Θ > 0, we solve the equation

( ) (1 ) 0 1 / gx G x x x x B αδ = -+ - = + Θ +
This leads to solving a third-order polynomial 0 0 1 2 2

3 3 = + + + a x a x a x a where Θ + Θ = - + Θ + Θ + = - Θ + + = = αδ αδ 0 1 2 3 / 1 / / / 1 1 / 1 a g B a B g B B a B a Since G is always concave, 2 3 3 2 2 / ''( ) 0 ( ) (1 / ) g B G x x x B αδ Θ = - - < + Θ +
model (3a) can have at most two different positive equilibria. We used numerical simulations in MATLAB to find out that there is either no positive equilibrium or that they are just two of which the lower is unstable (Allee threshold) and the upper is locally stable (carrying capacity).

A c c e p t e d m a n u s c r i p t 21

The case Θ = 0 (no mate-finding Allee effect, only Allee effect due to predation) allows for more detailed analysis. Straightforward yet tedious algebra gives the results summarized in Table A1; the areas I-IV are shown in Fig. A1.

For model (3b) and Θ > 0, the analysis follows the same lines as above with analogous results. The results for Θ = 0 are summarized in Table A2; the areas I-V are shown in Fig. A2.

For both models (3a) and (3b) and Θ > 0, the extinction equilibrium x 0 can also be unique and globally stable if the mate-finding Allee effect and/or the Allee effect due to predation are sufficiently strong; for α = 0 or B 0 or λ = 0, this happens when Θ > Figure 5: Parameter regions yielding the five different scenarios for populations with predation limit (A) and the three scenarios for populations without predation limit (B).

Θ lim = ( ) 2 1 g - .
See also Table 3 and the main text for details. A1.

Figure A2: Numbers and stability of equilibria of model (3b) in the λ-H parameter space for Θ = 0. Areas I-V explained in Table A2.

  respond numerically to changes in abundance of the focal population; b is the per capita birth rate, d is the intrinsic mortality at low densities, K scales the environmental carrying capacity, M(N) represents a mate-finding Allee effect or reduction of the reproductive rate bN due to difficulties in finding mates at low (N) is a predator functional response, and P represents a constant predator population density. The parameters are summarized together with other symbols in Table

  no mate-finding Allee effect, and predation ceases (α = 0, B = 0, and λ = 0) for max Θ ≥ Θ when reproductive activities become so suppressed that predators can no longer exploit the prey. The parameter z quantifies the overall predation risk associated with reproduction. Values of z < 1 (convex trade-offs) correspond to 'cheap' reproductive activities, for which the predation risk remains low until Θ is relatively small, i.e. until the mate-finding Allee effect is relatively weak due to efficient prey reproductive activities. On the other hand, values of z > 1 (concave trade-offs) correspond to 'costly' reproduction, for which the predation risk remains wide range of effort associated with reproductive activities, e.g. because courting or gravid individuals are more conspicuous (Fig.1).

  if Θ > Θ lim , while populations characterized by trade-off curve C go extinct for intermediate values of Θ between 0 and Θ max (curve C in Fig.

  Θ are typically 'the safest' and the population goes extinct for high, low, and both high and low levels of reproductive activity in scenarios C, D and

  taxa subject to the former include mostly vertebrates such as fishes (salmonids and others) and mammals (elk, caribou or pronghorn) but only a few invertebrates, including the crown-of-thorns starfish Acanthaster plancii or the rock lobster Jasus edwardsii. On the contrary, Allee effects due to mate limitation, including the mate-finding Allee effect, seem to occur mostly in invertebrates with the examples spanning molluscs (queen conch, abalone) and insects (gypsy moth, Glanville fritillary); only a few studies have identified this effect in fishes, mammals and birds.

  Haliotis sorenseni, in which individuals have reduced fertilization efficiency at small density and where some populations suffer overexploitation equivalent to an Allee effect due to predation. Additional research should reveal how frequent are such examples in other animals.

Figure 1 :

 1 Figure legends

Figure 2 :

 2 Figure 2: Equilibria of model (3a). Stable (solid lines) and unstable (dashed lines) equilibria of model (3a) as they depend on (A) the mate-finding Allee effect strength Θ and the predator attack rate α for g = 2, B = 0.01 and δ = 10 and on (B) the matefinding Allee effect strength Θ and the prey vulnerability B for g = 2, α = 0.05, δ = 10.Lines of equilibria in panel A range from α = 0 until 2.4 (step size 0.2; no equilibria beyond α = 2.55), and lines of equilibria in panel B start with B = 0.0001 and have step size 0.1 (lines for large B approach the leftmost line corresponding to B ∞).

Figure 3 :

 3 Figure 3: Population resilience (panels A and C) and relative population resilience for prey with inflexible reproductive behaviour (panels B and D) for model (3a), plotted as a function of the mate-finding Allee effect strength Θ and the predator attack rate α (panels A and B) and the mate-finding Allee effect strength Θ and the prey vulnerability B (panels C and D). Parameter values: A and B. g = 2, B = 0.01 and δ =10, so that Θ lim = 0.17 and α lim = 2.55. C and D. g = 2, α = 0.05 and δ = 10, so that Θ lim = 0.17 and there is no B lim . The thin lines are isolines of the (relative) population resilience, increasing from 0 to g-1 with stepsize 0.05 (indicated by thick arrow).These two pairs of panels (A with B and C with D) represent two generic outcomes that all models analyzed in this article produce.

Figure 4 :

 4 Figure 4: Population resilience under reproduction-predation risk trade-off. Panel A: The population resilience R as a function of Θ and α with examples of three trade-off curves marked A, B and C. Thin black lines = isolines of population resilience, increasing from 0 to g-1 with stepsize 0.05 (indicated by thick arrow). Populations cannot persist above the viability limit = the isoline connecting the points [Θ lim , 0] and [0, α lim ] at which R = 0; R is maximized for Θ = 0 and α = 0 where R = g-1. Panel B: The population resilience R as a function of Θ along the three trade-off curves shown in panel A. Trade-offs A and B yield the same resilience for Θ > Θ max = 0.08 and trade-offs A, B and C for Θ > Θ max = 0.12. See the main text for details. Parameter

  = 2, B = 0.01, δ = 10; trade-off curves: Θ max = 0.08, α max = 1, z = 0.4 (A), Θ max = 0.08, α max = 1, z = 3 (B), and Θ max = 0.12, α max = 1.8, z = 3 (C).

Figure 6 :

 6 Figure 6: Examples of the five scenarios of the impact of reproduction-predation risk trade-off between the mate-finding Allee effect strength Θ and the attack rate α on the population resilience R for the AR model. Panels A-E represent examples from the equally marked areas in Fig. 5A. Parameter values: g = 2, B = 0.01 and δ = 10, so thatΘ lim = 0.17 and α lim = 2.55; A: Θ max = 0.08 < Θ lim and α max = 1 < α lim ; B: Θ max = 0.08 < Θ lim and α max = 1.8 < α lim ; C: Θ max = 0.18 > Θ lim and α max = 1 < α lim ; D: Θ max = 0.08 < Θ lim and α max = 2.6 > α lim ; E: Θ max = 0.18 > Θ lim and α max = 2.6 > α lim .

Figure 7 :

 7 Figure 7: Examples of the five scenarios of the impact of reproduction-predation risk trade-off between the mate-finding Allee effect strength Θ and the attack rate α on the relative population resilience ρ i for prey with inflexible reproductive behaviour for the AR model. Legend as in Fig. 6.

Figure A1 :

 A1 Figure A1: Numbers and stability of equilibria of model (3a) in the α-B parameter space for Θ = 0. Areas I-IV explained in TableA1.

  Figure 1

Acknowledgements

We thank two anonymous reviewers for valuable suggestions on an earlier version of the manuscript. We acknowledge funding by the Institute of Entomology (Z50070508), Grant Agency of the Academy of Sciences of the Czech Republic (KJB600070602) and Faculty of Science (MSM6007665801).


Loris Nanni
email: loris.nanni@unibo.it

Jian-Yu Shi
email: jianyushi@nwpu.edu.cn

Sheryl Brahnam
email: sbrahnam@missouristate.edu

Alessandra Lumini
email: alessandra.lumini@unibo.it

Protein classification using texture descriptors extracted from the protein backbone image

Keywords: protein classification, texture descriptors, primary structure, local binary patterns, Radon transform, Haralick features, support vector machines

In this work we propose a method for protein classification that combines different texture descriptors extracted from the 2-D distance matrix obtained from the 3-D tertiary structure of a given protein. Instead of considering all atoms in the protein, the distance matrix is calculated by considering only those atoms that belong to the protein backbone. The positive results reported in this paper offer further experimental confirmation that the distance matrix contains sufficient information for describing a protein. Moreover, we show that combining features extracted from the primary structure with features extracted from the distance matrix increases the performance of our classification system. We demonstrate this finding by comparing the performance of an ensemble of classifiers that uses the combined features. The classifiers used in our experiments are support vector machines and random subspace of support vector machines. The experimental results, validated using three different datasets (protein fold recognition, DNA-binding proteins recognition, biological processes and molecular functions recognition) along with different texture feature extraction methods (variants of local binary patterns, radon feature transform based approaches, and Haralick descriptors) demonstrate the effectiveness of the proposed approach. Particularly

Introduction

An important area of research involves finding good methods for extracting a set of features from a protein [START_REF] Chou | Review: Prediction of protein structural classes[END_REF]. There are two general views on how extraction should be accomplished. Based on the wide assumption that structural features are closely related to sequence composition (Krissinel, 2007;Bastolla et al., 2008), one popular approach, called the indirect representation of protein spatial structure, extracts features from a sequence to perform classification. Indirect representation can be organized mainly into two types: one based on the statistical analysis of amino acid residues [START_REF] Shi | Using Pseudo Amino Acid Composition to Predict Protein Subcellular Location: Approached with Amino Acid Composition Distribution[END_REF] and the other based on amino acid indices (Cai et al., 2002;[START_REF] Shi | Using Pseudo Amino Acid Composition to Predict Protein Subcellular Location: Approached with Amino Acid Composition Distribution[END_REF]. Probably the most famous method for extracting a set of features from the amino-acid sequence is Chou's pseudo amino acid (PseAA) composition approach, and several variants of this method have been proposed in the literature, including hydropathy scales [START_REF] Chou | Using amphiphilic pseudo amino acid composition to predict enzyme subfamily classes[END_REF], physicochemical distance [START_REF] Chou | Prediction of protein subcellular locations by incorporating quasi-sequence-order effect[END_REF], digital code [START_REF] Gao | Using pseudo amino acid composition to predict protein subcellular location: approached with Lyapunov index, Bessel function, and Chebyshev filter[END_REF], complexity factor [START_REF] Xiao | Using complexity measure factor to predict protein subcellular location[END_REF]Xiao et al., 2006b), digital signal [START_REF] Chou | Digital coding of amino acids based on hydrophobic index[END_REF], Fourier low-frequency spectrum (Liu et al., 2005), cellular automata (Xiao et al., 2006a), and genetic programming (Nanni & Lumini, 2008b). As summarized in recent comprehensive reviews, (Chou & Shen, 2007;[START_REF] Chou | Pseudo amino acid composition and its applications in bioinformatics, proteomics and system biology[END_REF], since the concept of pseudo amino acid composition (PseAAC) was proposed [START_REF] Chou | Prediction of protein cellular attributes using pseudo amino acid composition[END_REF], it has provided a very flexible mathematical frame for investigators to incorporate their desired information into the representation of protein samples. According to its original definition, the PseAAC is actually formulated by a set of discrete numbers [START_REF] Chou | Prediction of protein cellular attributes using pseudo amino acid composition[END_REF] as

A c c e p t e d m a n u s c r i p t

3 long as it is different from the classical amino acid composition (AAC) and it is derived from a protein sequence that is able to harbor some sort of its sequence order and pattern information, or able to reflect some physicochemical and biochemical properties of the constituent amino acids.

The PseAAC approach has also been widely used to deal with many protein-related problems and sequence-related systems (see, e.g., (Chen et al., 2009;[START_REF] Ding | Prediction of cell wall lytic enzymes using Chou's amphiphilic pseudo amino acid composition[END_REF]Ding and Zhang, 2008;[START_REF] Fang | Predicting DNA-binding proteins: approached from Chou's pseudo amino acid composition and other specific sequence features[END_REF][START_REF] Georgiou | Use of fuzzy clustering technique and matrices to classify amino acids and its impact to Chou's pseudo amino acid composition[END_REF][START_REF] Gonzalez-Diaz | Proteomics, networks, and connectivity indices[END_REF]Jiang et al., 2008;[START_REF] Fang | Predicting DNA-binding proteins: approached from Chou's pseudo amino acid composition and other specific sequence features[END_REF]Lin, 2008;Lin et al., 2008;[START_REF] Ding | Prediction of cell wall lytic enzymes using Chou's amphiphilic pseudo amino acid composition[END_REF]Qiu et al., 2009;Zeng et al., 2009;Zhang and Fang, 2008;[START_REF] Shi | Using Pseudo Amino Acid Composition to Predict Protein Subcellular Location: Approached with Amino Acid Composition Distribution[END_REF]Zhou et al., 2007) and a long list of PseAAC-related references cited in a recent review [START_REF] Chou | Pseudo amino acid composition and its applications in bioinformatics, proteomics and system biology[END_REF]). As summarized in [START_REF] Chou | Pseudo amino acid composition and its applications in bioinformatics, proteomics and system biology[END_REF], until now 16 different

PseAAC modes have been used to represent the samples of proteins for predicting their attributes.

Each of these modes has its own advantage and disadvantage. Other new developments for predicting various protein features based on sequence information are: [START_REF] Ding | Prediction of cell wall lytic enzymes using Chou's amphiphilic pseudo amino acid composition[END_REF]Xiao et al., 2008a;Xiao et al., 2008b;Xiao et al., 2009a;Xiao et al., 2009b;Xiao et al. 2009c).

In contrast to the indirect approach is the view that features should be extracted directly from an analysis of the protein's spatial structure. This direct approach to feature representation can be grouped into three general types: one based on the spatial atom distribution [START_REF] Daras | Three-Dimensional Shape-Structure Comparison Method for Protein Classification[END_REF], a second on its topological structure [START_REF] Lumini | Voronoi and Voronoi-Related Tessellations in Studies of Protein Structure and Interaction[END_REF], and a third on its geometrical shape (Sayre & Singh, 2008).

Generally, the indirect representation is lower in computational cost but provides a higher dimensional feature set whereas the direct representation is higher in computational cost but provides a lower dimensional feature set. While the lower computational cost involved in the indirect approach is desirable, the higher dimensional representation requires the application of the most advanced techniques in pattern recognition, for example, building ensembles of classifiers for improving the performance of stand-alone methods [START_REF] Chou | Predicting protein-protein interactions from sequences in a hybridization space[END_REF]Sarda et al, 2005;Nanni & Lumini, 2006;Nanni & Lumini, 2008). (LTP) [START_REF] Tan | Enhanced local texture feature sets for face recognition under difficult lighting conditions[END_REF]. With DLP, the most frequent rotation invariant patterns are selected.

LTP is proposed for obtaining a noise robust texture descriptor. It is based on encoding the gray level difference d between a pixel x and its neighborhood u by 3 values.

A further improvement of performance is obtained by employing a "supervised" random subspace of classifiers where each bin of the histogram has a probability of belonging to a given subspace according to its occurrence frequencies in the training data (Nanni et al., 2009). Another approach for improving classification performance is to extract the texture descriptors not only from the whole distance matrix but also from some selected sub-windows.

In addition to investigating the performance of various texture descriptors, we build on the idea that descriptors based on different extraction notions give complementary information. We do this by combining a direct (Chou's amino acid) descriptor with an indirect representation (protein spatial structure features extracted from the distance matrix) using the sum rule. Our investigation shows that the classifiers trained with features extracted from the amino-acid sequence are partially independent to the classifiers trained with features extracted from the distance matrix. The experimental results show that the proposed ensemble of classifiers combining direct and indirect descriptors outperforms stand-alone approaches.

The fusion approaches are a well known technique for improving the performance of the standalone methods, remarkable successes of using fusion approaches are: (Chou and Shen, 2008b) protein subcellular location prediction, (Chou and Shen, 2008a) protease type prediction, (Chou and , 2009b;Shen et al., 2009) protein folding rate prediction, (Chou and Shen, 2007b) signal peptide prediction, (Chou and Shen, 2007a) membrane protein type prediction.

The remainder of this paper is organized as follows. In section 2, we introduce our proposed approach and the feature extraction methods investigated in this paper. In section 3, we report experimental results obtained on three benchmark databases. Finally, in section 4, we summarize results and draw a few conclusions.

Proposed approach

One objective in this work is to explore new methods for extracting features from the distance matrix that works well on several benchmark datasets. In addition, we want to investigate fusion using a standard method for extracting features from the primary sequence of the protein. The protein descriptor used in our experiments is Chou's well-known pseudo amino acid descriptor [START_REF] Chou | Using amphiphilic pseudo amino acid composition to predict enzyme subfamily classes[END_REF].

In [START_REF] Shi | Using Texture Descriptor and Radon Transform to Characterize Protein Structure and Build Fast Fold Recognition[END_REF], the authors show that Haralick features and the Radon transform produce a good texture descriptor for the distance matrix. In this paper, we compare several variants of the LBP, which we use as a new texture descriptor. We demonstrate that it is well suited for the distance matrix.

To improve performance further, we use a random subspace of support vector machines as the classifier, and we extract the features not only from the whole distance matrix but also from the II, III, and IV quadrant of the whole image (see figure 1). For each image a different classifier is trained and the results are combined using the sum rule. The architecture of our best performing system (random subspace of support vector machines) is presented in figure 2.

A general description of each step in our classification experiments is provided below.

Extracting features from the distance matrix of the protein backbone

Protein is a composite of 20 types of amino acid residues. Various physicochemical properties and different counts and sequenced orders of these residues are the keys for deciding and producing the diversity of protein spatial structures. Unfortunately, how these keys work together is not fully understood. This fact brings out the difficulty of describing, analyzing, and characterizing protein structure. Instead of considering all atoms, many researchers use the C α atoms of protein to characterize the whole protein structure. This protein backbone reflects the topology and the folding of protein [START_REF] Taylor | Protein structure alignment[END_REF]). An effective representation of the backbone is the distance matrix (DM). It contains sufficient information of the proteins structure as the original 3D backbone structure can be reconstructed from DM using distance geometry methods [START_REF] Timothy | The theory and practice of distance geometry[END_REF].

Given a protein P i as first step we need to extract its backbone: it is described as a vector

{ }

,1 ,2 , , , ,

i i i i N B α α α = Coor Coor Coor
, where , i n α

Coor is coordinates vector of the n th C α atom. The distance matrix is defined as the matrix

( ) ( ) { } , ,
, ,

i i i p q DM dm p q dist α α = = Coor Coor where ( ) dist ⋅ is
simply the Euclidean distance between the two set of coordinates (considered as a vector) and

1 , p q N ≤ ≤ . 1
Since DM maintains sufficient 3-D structural information, similar protein backbones are expected to have such distance matrices with similar properties. In our model, DM is regarded as a grayscale image. It is interesting to note that features extracted from DM are invariant to rotation and translation. 1 The matlab code for extracting the distance matrix is available at http://bias.csr.unibo.it/nanni/DM.zip In (Liao et al., 2009), better performance (compared with using the uniform patterns) is obtained when the patterns that represent K% (K=90% in their work) of the whole pattern occurrences in the training data are selected. To make LBP more robust to noise, [START_REF] Tan | Enhanced local texture feature sets for face recognition under difficult lighting conditions[END_REF] used Local Ternary Patterns (LTP). In LTP the difference d between x and its neighborhood u is encoded by 3 values according to a threshold τ (here τ=3):

1 if u ≥ x + τ ; -1 if u ≤ x -τ ; else 0.
The ternary pattern is then split into two binary patterns by considering its positive and negative components. The histograms that are computed from the binary patterns which are concatenated to form the feature vector.

In our experiments, we obtain the best performance by combining the idea of dominant LBP with LTP. For this purpose, we modify the original LBP code found at http://www. ee.oulu.fi/mvg/page/lbp_matlab2 .

A c c e p t e d m a n u s c r i p t

9 Furthermore, when we use dominant LTP as feature extractor, we do not use the standard random subspace method but rather a "supervised" approach. The standard random subspace method [START_REF] Ho | The random subspace method for constructing decision forests[END_REF] modifies the training dataset by generating K (where K=100 in this paper) new training sets containing only a random subset of 60% of the features. Classifiers are then trained on these modified training sets. In "supervised" random subspace, each feature is not randomly selected but has a probability of being chosen that is based on the occurrence frequencies in the training data. Given x i, the sum of the occurrence frequency of the i-th bin of the histogram in the training data, the probability of the i-th bin to be chosen is given by x i /∑ i x i. In "supervised" random subspace, the scores of the set of classifiers are combined using the sum rule.

It should be noted that the "supervised" random subspace method is only used with the dominant LBP/LTP features. When random subspace is coupled with other texture descriptors, the standard random subspace is used.

Radon transform

This feature extraction method is performed by selecting the first discrete cosine transform The Radon transform is the projection of the image intensity along a radial line oriented at a specific angle (in this paper we use an angle of 45°). We have also tested other methods for selecting discrete cosine coefficients 3 (DCT), e.g., the coefficients with higher variance in the training set and the "discriminative power" method proposed in [START_REF] Dabbaghchian | Feature ex-traction using discrete cosine transform and discrimination power analysis with a face recognition technology[END_REF]. But in our tests, the best performance is obtained by selecting the first discrete cosine coefficients. In figure 3 an example of this feature extraction is reported.

3 For both Radon and DCT we have used the official matlab function (radon.m and dct.m) The Haralick texture features descriptor was proposed 30 years ago by [START_REF] Haralick | Statistical and structural approaches to texture[END_REF]. It is based on the spatial gray level dependence matrices (SGLD), or co-occurrence matrix. Given an image with N gray levels, the SGLD matrix at angle θ is a matrix of size N×N. Each element in the matrix is a count of the total number of pairs of gray levels i and j at a distance d along the direction θ.

Thirteen features are calculated4 from a SGLD matrix at a fixed angle θ: energy, correlation, inertia, entropy, inverse difference moment, sum average, sum variance, sum entropy, difference average, difference variance, difference entropy, and two information measures of correlation.

In this work we test two different feature sets extracted using the Haralick's method:

• HARA1, concatenation of the Haralick features extracted by considering two angles ( 0° and 90°), with d=1;

• HARA3, concatenation of the features extracted by considering four angles (0°, 45°, Chou's pseudo amino acid feature vector a different support vector machine is trained. These 50 classifiers are then combined using the sum rule.

Classification system

For the classifier, we have used the well-known support vector machine (SVM) [START_REF] Cristianini | An introduction to Support vector machines and other kernelbased learning methods[END_REF]. It is a two class classifier which aims at finding the hyperplane that separates the training patterns of two classes by maximizing the distance between the hyperplane and the two classes. When it is not possible to find a linear decision boundary, a kernel function can be used to project the data onto a higher-dimensional feature space where a hyperplane separating the two classes can be found. Some typical kernels used in SVM include polynomial kernels and radial

A c c e p t e d m a n u s c r i p t

13 SVM-based machine learning algorithm was used in predicting protein subcellular location [START_REF] Chou | Using functional domain composition and support vector machines for prediction of protein subcellular location[END_REF], membrane protein type (Cai et al., 2003a;Cai et al., 2004a), protein structural class (Cai et al., 2002a), specificity of GalNAc-transferase (Cai et al., 2002b), HIV protease cleavage sites in protein (Cai et al., 2002c), beta-turn types (Cai et al., 2002d), protein signal sequences and their cleavage sites (Cai et al., 2003b), alpha-turn types (Cai et al., 2003c), catalytic triads of serine hydrolases (Cai et al., 2004b), B-cell epitope prediction (Chen et al., 2007).

Datasets

All experiments were performed using the following datasets: Protein fold recognition, DNAbinding proteins, and the GO dataset. Each of these datasets and the testing protocols are briefly described in this section.

Protein fold recognition (FOLD)

The fold database used in our experiments is derived from the work of [START_REF] Ding | Multi-class protein fold recognition using support vector machines and neural networks[END_REF]. It contains a training set and a testing set that contain 313 and 385 proteins. The sequence similarities are less 35% and 40% respectively, and the class numbers are both 27. The training set is used to build the classifier models, and we independently used the testing set to evaluate performance as this testing protocol is widely used in the literature for this dataset. The whole database can be downloaded from http://ranger.uta.edu/~chqding/protein/. Some sample of distance matrices extracted from the proteins of this dataset are reported in figure 4.

DNA-binding proteins (DNA)

This is the dataset is reported in [START_REF] Fang | Predicting DNA-binding proteins: approached from Chou's pseudo amino acid composition and other specific sequence features[END_REF] and contains 118 DNA-binding Proteins and 231 Non-DNA-binding proteins. These proteins have less than 35% sequence identity between each pair. DNA-binding proteins are proteins that are composed of DNA-binding domains and thus have a specific or general affinity for either single or double stranded DNA. Sequence-specific DNA-binding proteins generally interact with the major groove of B-DNA.

Some sample of distance matrices extracted from the proteins of this dataset are reported in figure 5. For this database we used the ten-fold cross validation protocol.

GO dataset (GO)

This dataset was reported in (Nanni et al., 2009b). It was created by collecting proteins according to GO annotations, distinguishing between the biological processes "immune response" (33 proteins), "DNA repair" (43 proteins), and between the molecular functions "substrate specific transporter activity" (39 proteins) and "signal transducer activity" (53 proteins). The presence of highly similar proteins in the same class was avoided by removing sequences having more than 30% identity. We randomly extracted 20% of the proteins for building the testing set, and this procedure was repeated 50 times. The results were then are averaged.

Some sample of distance matrices extracted from the proteins of this dataset are reported in figure 6.

Experimental results

Performance is evaluated using the area under the Receiver Operating Characteristic (ROC)

curve. The area under the ROC curve (AUC 6 ) [START_REF] Fawcett | ROC Graphs: Notes and Practical Considerations for Researchers[END_REF]) is a scalar measure to evaluate performance which can be interpreted as the probability that the classifier will assign a higher score to a randomly picked positive sample than to a randomly picked negative sample. In the GO dataset, which is a four class problem, the AUC is calculated using the one versus all approach (a given class is considered as "positive" and all the other classes are considered as "negative") and the average AUC is reported.

In the following tests, Table 1, the AUC is reported as obtained using the following methods:

• [START_REF] Shi | Using Texture Descriptor and Radon Transform to Characterize Protein Structure and Build Fast Fold Recognition[END_REF], the method reported in [START_REF] Shi | Using Texture Descriptor and Radon Transform to Characterize Protein Structure and Build Fast Fold Recognition[END_REF], where few selected Radon features and Haralick feature are extracted from the distance matrix;

• PSEAA, the Chou's pseudo amino acid composition method explained in section 2.5;

• LBP, standard LBP with P=16 and R=2, the entire DM is used;

• LTP, standard LTP with P=16 and R=2, the entire DM is used;

• DLBP, dominant LBP with P=16, R=2 and K=90%, the entire DM is used;

6 EUC is implemented as in dd_tools 0.95 davidt@ph.tn.tudelft.nl • RS-DLTP, a random subspace of DLTP, the entire DM is used;

• MI DLTP, fusion of four DLTPs, one extracted from the entire DM and the others extracted from three sub-windows of the DM;

• MI(RS-DLTP), fusion of four RS-DLTPs, one extracted from the entire DM and the others extracted from three sub-windows of the DM;

• RADON+DCT, the descriptor described in section 3.3, the entire DM is used and the first 20 DCT coefficients are retained;

• RS-RADON+DCT, the random subspace version of RADON+DCT;

• MI(RS-RADON+DCT), fusion of four RS-RADON+DCT, one extracted from the entire DM and the others from three sub-windows of the DM. The first 20 DCT coefficients are retained when the DM is used but with the sub-windows the first 10 coefficients are retained;

• HARA1, the descriptor described in section 3.4, only the entire DM is used;

• HARA3, the descriptor described in section 3.4, only the entire DM is used;

• MI(HARA3), fusion of four HARA3, one extracted from the entire DM and the others from three sub-windows of the DM.

• In each cell of table 1 when the texture descriptors are used, there are two values: the first is the AUC obtained when the original distance matrix is used, and second, a value between parenthesis, is the AUC obtained when the matrix is resized to 100×100 before the feature extraction step. Also, for MI(RS-DLTP) and MI HARA3, we use the images resized to 100×100 before the feature extraction step, while for MI(RS-RADON+DCT) we use the original images. From the results reported in table 1, the following conclusions can be drawn:

A c c e p t e

• Using the three sub-windows of the DM improves performance;

• The dominant LTP with the LBP/LTP texture descriptor works well in this particular application--the standard LBP/LTP obtains the worse performance;

• The features used in [START_REF] Shi | Using Texture Descriptor and Radon Transform to Characterize Protein Structure and Build Fast Fold Recognition[END_REF], a set of selected Haralick features and moments extracted from the Radon coefficients, works particularly well in the BINDING dataset and well enough in the FOLD dataset. In these datasets the structural configuration of the proteins is very important. We want to stress the low performance of PSEAA in the FOLD dataset. The GO dataset the features of [START_REF] Shi | Using Texture Descriptor and Radon Transform to Characterize Protein Structure and Build Fast Fold Recognition[END_REF] works poorly probably because of the low performance of the Radon based features in this difficult dataset.

A c c e p t e d m a n u s c r i p t 19

Notice that in our work we do not perform any feature selection. We prefer to use longer feature vectors since the feature selection could be dataset dependent. In our work we want to propose a general method that works well on several datasets. If a training set is used to select a set of features for a particular dataset, the performance would likely be improved.

In table 2, we report some tests performed for optimizing the parameters of the texture descriptors, as follows:

• Different values for the parameters τ and K of dominant LTP;

• MULTIRES, a combination between a dominant LTP with P=16 and R=2 and P=8 and R=1, both with τ=3;

• R X-Y, is the MI(RS-RADON+DCT) method where the first X coefficients are extracted from the entire DM, and the first Y coefficients are extracted from the sub-windows of DM.

• MI HARA3 d=X, is the method MI HARA3 where we change the value of the parameter d. Finally, we report the results of the following fusions by sum rule (notice that before the fusion the scores of each classifier are normalized to a mean of 0 and a standard deviation of 1):

A c c

• FUS1, fusion by sum rule among MI(RS-DLTP), MI(RS-RADON+DCT) and MI HARA3; worst in the structural classification problem (where the features extracted from the distance matrix work very well). Obviously, in this particular problem, the structural classification using the distance matrix brings more information than the amino-acid sequence.

FUS3 obtains the best performance in the FOLD dataset, while in the BINDING/GO dataset the best performance is obtained by FUS4 (where PSEAA has an higher weight).

The results also suggest that some ad-hoc fusion methods should be studied independently in each dataset. For example, in the BINDING dataset fusion by sum rule between MI HARA3 and PSEAA (the two best descriptors in that dataset) obtains an AUC of 92.52. If in the BINDING dataset, we combine using sum rule MI HARA3, PSEAA and [START_REF] Shi | Using Texture Descriptor and Radon Transform to Characterize Protein Structure and Build Fast Fold Recognition[END_REF], we obtain an an AUC of 93.58.

In order to validate the effectiveness of the presented method, we compare it with several methods in different literatures all of which used the same benchmark data sets and the same testing protocol of original paper [START_REF] Ding | Multi-class protein fold recognition using support vector machines and neural networks[END_REF], it is the dataset named FOLD in this work. The comparison is listed in Table 4. [START_REF] Ding | Multi-class protein fold recognition using support vector machines and neural networks[END_REF] proposed six kinds of features denoted by Hydrophobicity and Predicted secondary structure respectively. [START_REF] Chinnasamy | Protein structure and fold prediction using tree-augmented naive bayesian classifier[END_REF]Shi et al., 2006) are based on the same features CSHPVZ but different classifier systems. [START_REF] Huang | Hierarchical learning architecture with automatic feature selection for multiclass protein fold classification[END_REF] combines CSHPVZ with bigram-coded feature (B) and spaced bigram-coded feature(SB), (Lin et al., 2007) does the same work as [START_REF] Huang | Hierarchical learning architecture with automatic feature selection for multiclass protein fold classification[END_REF] but improves the classifier system by the technique of data fusion.

Method

Accuracy (%) [START_REF] Ding | Multi-class protein fold recognition using support vector machines and neural networks[END_REF] 56.50 [START_REF] Chinnasamy | Protein structure and fold prediction using tree-augmented naive bayesian classifier[END_REF] 58.18 (Shi et al., 2006) 61.04 [START_REF] Huang | Hierarchical learning architecture with automatic feature selection for multiclass protein fold classification[END_REF]

Conclusion and Discussion

This paper focused on the study of texture descriptors for training an ensemble of machine learning algorithms for protein classification. The texture descriptors are extracted from the 2-D distance matrix obtained from the 3-D tertiary structure of a given protein.

Based on an analysis of prior research, we propose a new method based on a fusion of three texture descriptors and on the pseudo Chou's amino acid descriptor. Moreover, to improve the The ensemble proposed in this work has been tested on three datasets. The experimental results show that the proposed ensemble of classifiers outperforms stand-alone approaches.

Particularly interesting are the results in the classification of 27 types of structural properties.

The best practical finding revealed in this work is that texture descriptors extracted from the 2-D distance matrix and amino acid descriptors should be combined to obtain a very reliable method of classification.

As further future work we want to test more texture descriptors, some preliminary tests have shown interesting results. We have tested: an holistic method based on the neighborhood preserving embedding method (NPE) 7 (He et al., 2005), which is a subspace learning algorithm aimed at preserving the global Euclidean structure of the space, obtaining a 40% accuracy in the FOLD dataset; a recent Gabor based descriptors (Guo et al., 2009) obtaining a 66% accuracy in the FOLD dataset. Both these results are obtained considering also the sub-windows of the distance matrix and the random subspace as classifier.

Since user-friendly and publicly accessible web-servers represent the future direction for developing practically more useful predictors (Chou and Shen, 2009a), we shall make efforts in our future work to provide a web-server for the method presented in this paper.

  we apply new pattern recognition techniques to the indirect representation of protein features by examining texture descriptors extracted from the 2-D distance matrix of different protein classification datasets. We validate the use of texture descriptors by testing many variants of state-of-the-art texture descriptors. In particular, we examine variants of the well-known local binary patterns. Our experiments show that the best performance is obtained when the idea of dominant local binary patterns (DLB) (Liao et al., 2009) is combined with local ternary patterns

Figure 1 .

 1 Figure 1. Subwindows of the whole distance matrix.

Figure 2 .

 2 Figure 2. Proposed system for protein classification

  The basic idea of LBP is to examine the joint distribution of gray scale values of a circularly symmetric neighbor set of P pixels around a pixel x on a circle of radius R. The LBP histogram of dimension N (usually N=P+2) is obtained considering all the pixels of a given image. The difference d between x and its neighborhood u is encoded by 2 values: A pattern that contains at most two bitwise 0 to 1 or 1 to 0 transitions (circular binary code) is called uniform patterns (e.g., 11111111, 00000110 or 10000111 are all uniform patterns). In uniform LBP only these distributions are considered in the histogram (where a single bin contains all the non-uniform patterns).

(

  duda et al., 2001) coefficients from the Radon Transform of the image (Kourosh & Hamid, 2005).

Figure 3 .

 3 Figure 3. Radon+DCT feature extraction.

  et al., 2009) a sequence-based algorithm combining the augmented Chou's pseudo amino acid composition based on auto covariance is presented. A set of pseudo amino acid based features are extracted from a given protein as the concatenation of the 20 standard amino acid composition values and m (where m=20) values, reflecting the effect of sequence order (where m is a parameter denoting the maximum distance between two considered amino acids i,j): k) denotes the index of the amino acid in the k th position of the protein, Len is the length of the protein, d denotes the selected physicochemical property, and the function index(i,d) returns the value of the property d for the amino acid i. M d and V d are normalization factors denoting the average and the variance of the physicochemical property d on the 20 amino acids: We create 50 different Chou's pseudo amino acid feature vectors using 50 different psychochemical properties extracted from the AAindex (Kawashima & Kanehisa, 2000). For each

  basis function kernels. It should be noted that the features used for training SVM are linearly normalized to [0 1]. In our work we use the OSU matlab toolbox (http://sourceforge.net/projects/svm/).

Figure 4 .

 4 Figure 4. Examples of different classes of the FOLD dataset.

Figure 5 .

 5 Figure 5. Examples of different classes of the BINDING dataset.

Figure 6 .

 6 Figure 6. Examples of different classes of the GO dataset.

•

  DLTP, dominant LTP with P=16, R=2 and K=90%, the entire DM is used;

•

  FUS2, fusion by weighted sum rule among MI(RS-DLTP), MI(RS-RADON+DCT) and MI HARA3. The weights of MI(RS-DLTP) and MI(RS-RADON+DCT) are 0.25, while the weight of MI HARA3 is 1; • FUS3, fusion by sum rule among PSEAA, MI(RS-DLTP), MI(RS-RADON+DCT) and MI HARA3; • FUS4, fusion by weighted sum rule among PSEAA, MI(RS-DLTP), MI(RS-Radon+Dct) and MI HARA3. The weight of PSEAA is 5, while the weights of the other methods are 1.

  ,H,P,V and Z respectively. Letter C is just the popular amino acid composition, while the left five letters indicate the features of Polarity, Polarizability, Normalized Van Der Waals volume,

  extract the texture descriptors not only from the entire distance matrix but also from some selected sub-windows.

Table 1 .

 1 Comparison among different methods.

	d	m a n u s c r i p t

Table 2 .

 2 Parameters optimization.

			DATASETS
		FOLD	BINDING	GO
	τ=3 K=90%	75.58	82.00	62.20
	τ=1 K=90%	70.17	79.67	59.24
	τ=5 K=90%	74.03	81.56	60.37
	τ=3 K=80%	72.99	77.80	60.52
	τ=3 K=85% τ=3 K=95% MULTIRES R 25-10 R 40-20 R 15-5 MI HARA3 d=1 MI HARA3 d=2 MI HARA3 d=3 From the results reported in table 2, it is clear that the best configurations are R 40-20 for 74.55 80.28 60.62 74.81 83.75 62.07 73.25 80.78 61.31 71.17 79.29 55.21 70.65 81.06 59.95 66.23 80.27 51.60 71.69 87.32 61.10 72.47 85.25 63.64 PARAMETERS 74.03 82.98 62.76 a n u s c r i p t m e p t e d RADON+DCT, d=1 for MI HARA3, and τ=3 and K=90% for dominant LTP.
				20

Table 3 .

 3 Fusion approaches.Looking at table 3, it is clear that the best approach is different for the FOLD dataset and for the BINDING/GO dataset. This is due to the different performance of PSEAA. It works better than a texture descriptor in the general classification problems in the BINDING or GO dataset but works

	Method FOLD BINDING GO
	FUS1	80.78	87.97	64.44
	FUS2	76.62	89.41	63.59
	FUS3	81.04	91.06	69.16
	FUS4	64.68	92.16	71.53
	A c c e p t e d	m a n u s c r i p t
				21

Table 4 .

 4 Other comparisons in the FOLD dataset.The results reported in table 4 demonstrate that our proposed system (FUS3) outperforms other methods with the highest accuracy of classification.

	65.50

we have extracted the uniform patterns with getmapping011(16,'riu2') while for dominant LBP/LTP we use getmapping011(16,'u2')

Implemented as in Haralick Texture Features Matlab Toolbox v0.1b shalinig@ece.utexas.eduwww.bme.utexas.edu/reasearch/informatics

 A c c e p tXiao, X., Wang, P., and Chou, K.C., 2008b. Predicting protein structural classes with pseudo amino acid composition: an approach using geometric moments of cellular automaton image. J

Theor Biol 254, 691-696.