PAGES: 641-648 DOI: Full paper
Dynamics of Experimental Populations of Native and Introduced Blowflies (Diptera: Calliphoridae): Mathematical Modelling and the Transition from Asymptotic Equilibrium to Bounded Oscillations

WAC Godoy, CJ Von ZubenI +, SF dos ReisII + ++, FJ Von ZubenIII +

Departamento de Parasitologia, IB, Universidade Estadual Paulista, 18618-000 Botucatu, SP, Brasil
ICurso de Pós-Graduação em Ciências Biológicas, Universidade Estadual Paulista, 13506-900 Rio Claro, SP, Brasil
IIDepartamento de Parasitologia, IB
IIIDepartamento de Computação e Automação Industrial, FEE Universidade Estadual de Campinas, Caixa Postal 6109, 13083-970 Campinas, SP, Brasil


The equilibrium dynamics of native and introduced blowflies is modelled using a density-dependent model of population growth that takes into account important features of the life-history in these flies. A theoretical analysis indicates that the product of maximum fecundity and survival is the primary determinant of the dynamics. Cochliomyia macellaria, a blowfly native to the Americas and the introduced Chrysomya megacephala and Chrysomya putoria, differ in their dynamics in that the first species shows a damping oscillatory behavior leading to a one-point equilibrium, whereas in the last two species population numbers show a two-point limit cycle. Simulations showed that variation in fecundity has a marked effect on the dynamics and indicates the possibility of transitions from one-point equilibrium to bounded oscillations and aperiodic behavior. Variation in survival has much less influence on the dynamics.

Chrysomya megacephala (F.) is a blowfly that originally occurs in Australasia, the Oriental and Paleartic regions, South Africa and Afrotropical islands (Baumgartner & Greenberg 1984, Smith 1986), whereas C. putoria (Wied.) ranges in distribution from Tanzania to Congo in Africa (Zumpt 1965, Smith 1986). These flies were first detected in South America around 1975 (Guimarães et al. 1978), and have since then become established in the Americas (Guimarães et al. 1979, Baumgartner & Greenberg 1984). These blowfies have dispersed rapidly throughout the continent and C. megacephala already has reached North America and has been reported in California (Greenberg 1988, Wells 1991). Blowflies of the genus Chrysomya have considerable medical and veterinary importance because they are known to produce myiasis in humans and other animals (Zumpt 1965, Baumgartner & Greenberg 1984) and may serve as mechanical vectors of enteric pathogens and parasites (Furlanetto et al. 1984, Wells 1991). These flies are also considered of potential forensic importance since they can breed in carrion exposed to environmental conditions (Catts & Goff 1992, Von Zuben et al. 1996).

The invasion of these economically important blowflies has apparently caused the sudden decline in population numbers of an ecologically similar species, Cochliomyia macellaria, which is native to the Americas (Guimarães et al. 1979, Prado & Guimarães 1982, Greenberg & Szyska 1984). This case of biological invasion, like many others, has classical features such as the rapid spread of the invaders which, in turn, can bring about the decline of native species at the local and macrogeographic scales (Hengeveld 1989, Lodge 1993). Biological invasions represent the outcome of complex interactions at the ecologic and genetic levels (Pimentel 1993, Brown 1993), and the dynamic behavior of populations remains as an important component in the evaluation of fundamental aspects for the success of invasions, such as survival and extinction rates in populations (Hengeveld 1989).

Recently, we initiated a research programme aimed at understanding the dynamics of equilibrium of native and introduced blowflies. Our approach has used mathematical models based on non-linear finite difference equations that revealed that the native and introduced blowflies differ markedly in their equilibrium dynamics (Godoy et al. 1993, Von Zuben et al. 1993, Reis et al. 1996). Application of the mathematical model using parameters derived from experimental populations showed that the two introduced species will form stable oscillations with numbers fluctuating over a three to four fold range in successive generations (Godoy et al. 1993, Von Zuben et al. 1993), whereas in the native species, C. macellaria, the dynamics is characterized by damping oscillations in population size leading to one fixed point equilibrium (Reis et al. 1996).

These results were obtained using fixed parameter values estimated from statistical analysis of data for fecundity and survival from experimental populations. Nevertheless, it is well known that changes in parameter values that control the growth rate in populations can change the dynamics from one-point equilibria to bounded oscillations and continuous chaos (May 1976, May & Oster 1976, Murray 1991). We address ourselves here to the problem of how changes in parameter values in the mathematical model affect the equilibrium dynamics of introduced and native blowflies. In what follows we will (1) describe the rationale for the mathematical approach taken in our modelling effort, which is based on the life-history of blowflies, (2) describe the stability properties of the mathematical model, (3) describe the equilibrium dynamics of native and introduced species, and (4) evaluate the sensibility of parameters that govern the stability of the mathematical model and how this extrapolates into changes in the equilibrium dynamics of experimental populations ofC. macellariaC. megacephala and C. putoria.



Laboratory populations of C. macellariaC. megacephala and C. putoria were founded from specimens collected in the campus of the Universidade Estadual de Campinas, Campinas, State of São Paulo, Brazil. Adult flies were maintained in laboratory conditions in cages (30 x 30 x 48 cm) covered with nylon at 25 ± 1o C and were fed water and sugar ad libitum. Adult females were fed fresh beef liver to permit the complete development of the gonotrophic cycle. The experiments were performed using the generation F2, which is progeny of one generation that had its life cycle completed in the laboratory. Exploitative competition among larvae, which is known to occur under natural conditions (de Jong 1976, Lomnicki 1988), was established in the laboratory by setting up 10 larval densities ranging from 200 to 2,000 larvae per vial at intervals of 200 for each species. 
For densities 200-1,200 two replicates were used, whereas for densities 1,400-2,000 only one replicate was used. Fecundity was measured counting the number of eggs per female, expressed as average daily egg output, based on the length of the gonotrophic cycle at 25oC (Avancini & Prado 1986, Linhares 1988). Maximum sample size for estimation of fecundity was 30 females per vial. Sample sizes smaller than 30 in some vials were due to either low immature survival rates or incomplete ovarian development. Survival was estimated as the number of adults emerging from each vial. Exponential regressions for fecundity and survival as functions of increasing larval density were fitted to the data for the three species using the regression procedure of SAS (SAS Institute 1988). Details of the mathematical procedures and simulations are given in the appropriate sections below.



Mathematical modelling of blowfly population dynamics - In natural populations of blowflies, the adults disperse in the search for substrates to feed and lay eggs (Hanski 1987). The substrates are discrete and ephemeral and the developing immatures experience high levels of competition for limited resources (Atkinson & Shorrocks 1981, Godoy et al. 1996). The amount of food consumed by larvae will determine the size of adults which is correlated with female fecundity, that is, larger blowfly females lay more eggs than smaller ones (Kamal 1958, Goodbrod & Goff 1990). The life-history of blowflies is thus characterized by generations that occupy patchy and discrete resources, and fecundity and survival of adults are modulated by density conditions that affect the larval stage. Therefore, the modelling of the population dynamics has to consider discrete generations and density dependence at the immature stage with a delayed effect on the survival and fecundity of adults. Prout and McChesney (1985) have developed a model that takes into account all these features. This model is based on a finite-difference equation that models the density-dependent dynamics of immatures, eggs or larvae, in succeeding generations, nt + 1 and nt, as a function of the decrease in fecundity (F) and survival (S) with increasing immature (n) density. The discrete-time population dynamics is then written as

nt + 1 = ½ F(ntS(ntnt (1)

where fecundity and survival are decreasing functions of the number of immatures. The factor ½ indicates that only half of the population are adult females that contribute eggs to the next generation. The non-linearity of (1) is given by the product F(ntS(nt). Using exponential functions, therecursion equation describing changes in populations numbers reads

nt + 1 = ½ FS*e-(f+s)n t nt (2)

where F* and S* are regression intercepts describing theoretical values of maximum fecundity and survival, respectively, and f and s are regression coefficients which estimate the dependence of fecundity and survival on increasing levels of larval density. The rationale for using exponential functions in (2) is explained below.

The equilibrium dynamics of populations obeying equation (2) can be assessed by the eigenvalue l calculated at k, which is the number of immatures at equilibrium, that is, nt+1 = nt = k. For exponential functions these two quantities are written, respectively, as

l = 1- ½ kF*f efk S(k)- ½ kS*s-sk F(k) (3)


k =1n F*S*/ 2 /(f+s). (4)

Stability properties of the mathematical model - In this section we show that the stability properties of the Prout and McChesney's model hinges solely on the product of maximum fecundity and survival. Only the main results are given here, and the technical details of the stability analysis can be found elsewhere (Teixeira et al. 1996). Equation (2) can be written for notational convenience as

nt+1=ae-bnt nt , (5)

where a = ½ F*S* and b = f + s. A steady-state solution n (equilibrium point) associated with equation (5) is defined to be the value that satisfies the relation

nt+1 = nt = n , (6)

so that no change occurs in subsequent generations (Murray 1991). Equation (5) has equilibrium points given by

n = a e-bnn , (7)

which implies that

n = 0 or n = 1na1/b(8)

The derivative of equation (5) with respect to nt, evaluated at the non-zero equilibrium point n = 1na1/b , yields

dnt+1=1-1na . dnt (9)

This equation indicates that the parameter 
a = ½ F*S* is the primary determinant of the dynamics. Thus, the stability of the equilibrium point n = 1na1/b depends only on a, which must satisfy 
|1- ln a| < 1 (Murray 1991), leading to

1 < a < e2. (10)

The evolution of stable equilibrium points in equation (5) as a function of a is shown in Fig. 1. In the interval 1 < a < e2there is one fixed point. Continually increasing a beyond a = e2 gives rise to a hierarchy of bifurcating stable equilibrium points that evolve to the chaotic regime. The dynamics of equation (5) is similar to that given by the logistic and Ricker maps used to model density-dependent population dynamics (May 1976, May & Oster 1976), and also displays the universal behavior of period-doubling route to chaos (Fig. 1), first discovered by Feigenbaum (1983). Equilibrium dynamics of native and introduced blowflies - The equilibrium dynamics of experimental populations of C. macellariaC. megacephala and C. putoria may be inferred using equation (2) employing regression estimates of dependence of fecundity and survival on the increasing number of immatures. Parameter estimates using exponential regressions are given in the Table. Exponential regressions were chosen because they provide a better fit than linear and hyperbolic regressions under similar experimental conditions (Godoy et al. 1993, Von Zuben et al. 1993, Reis et al. 1996). The use of exponential functions is also supported by the findings of Rodriguez (1989), which indicate that the decrease in fecundity as a function of density of immatures can be viewed biologically as a Poisson process, which is described by an exponential function.

Using estimates for F*, S*, f, and s in the Table with equation (2) we can construct the recurrence relations for C. macellariaC. megacephala and C. putoria, respectively, as follows:

nt + 1 = ½ [(18.33 e-0.000477nt ) (0.788 e-0.000985nt)] nt,

nt + 1 = ½ [(23.49 e-0.000624nt) (0.916 e-0.000148nt)] nt,

nt + 1 = ½ [(19.32 e-0.000569nt) (0.970 e-0.000135nt)] nt. (11)

The qualitative dynamic behavior of experimental populations in the three species can be appreciated interating the set of equations (11) through time, which shows that the qualitative dynamics differs noticeably between the native and the invading blowflies (Fig. 2). Cochliomyia macellaria shows damping oscillatory behavior leading to a one-point equilibrium, whereas in the invading species, C. megacephala and C. putoria, population numbers show bounded oscillations characterizing a two-point limit cycle. This distinct qualitative behavior can also be inferred from the eigenvalues (l). Using equations (3) and (4) and the parameters estimated from the data (Table) one obtains values of l as -0.9776, -1.3761, and -1.2399 for C. macellariaC. megacephala, and C. putoria, respectively. Eigenvalues less than 1 in module indicate a one-point equilibrium, whereas values larger than 1 indicate higher order cycles (Murray 1991).

Transition from asymptotic stable equilibrium to bounded oscillations and aperiodic behavior - The dynamic behaviors observed for the native and invading blowflies were derived from the application of fixed parameter values to the non-linear difference equation (2). The effect of parameter variation on the population dynamics of C. macellariaC. megacephala and C. putoria was investigated running simulations varying maximum fecundity and survival, F* and S*, respectively. All simulations were carried out with Matlab (Moler et al. 1987). In the simulations reported here, fecundity was arbitrarily allowed to vary up to a mean daily egg output of 40 eggs.

Bifurcation diagrams show that increasing values for fecundity do produce qualitative changes in the dynamics of the three species (Fig. 3). In the native blowfly, C. macellaria, the dynamics shifts from one-point equilibrium to bounded oscillations in population size, including two- and four-point limit cycles (Fig. 3). The two invading species, C. megacephala and C. putoria, whose dynamics was already oscillatory show increasing number of bounded cycles with increasing fecundity, and, for values of fecundity larger than 34 eggs, there is apparently a transition in the dynamics from bounded oscillations in population number to a region of aperiodic oscillations (Fig. 3). We point out that the upper limit of 40 eggs used in our simulations does not match exactly the values observed for these and other calliphorid species, since the maximum number of eggs can be either smaller or larger depending upon the species and the level of competition experienced by the larvae (Wijesundara 1957, Ribeiro 1992, Godoy et al. 1993, Von Zuben et al. 1993, Reis et al. 1994).

In the simulations involving bifurcations as a function of the variation of survival, S, the upper limit considered was 1.0, which represents the maximum viability of 100%. The effect of survival on the dynamics of blowflies is much less noticeable than that of fecundity, since no more than one bifurcation is realized in the three species (Fig. 4). Increasing values of survival do not change the qualitative dynamics of C. megacephala or C. putoria. On the other hand, values of survival larger than 0.8 place the native species, C. macellaria, in the two-point limit cycle zone.

Our analytical results indicated that the product of maximum fecundity and survival, FS*, is the primary determinant of the dynamics of populations obeying equation (2). Of these two parameters, fecundity was shown in the simulations to be the most important demographic parameter to bring about shifts in the dynamic behavior of experimental populations of native and introduced blowflies. Increasing values for fecundity place the native species, C. macellaria, in the region of bounded oscillations in population size. Higher order cycles and even aperiodic oscillations may arise in the introduced blowflies, C. megacephala and C. putoria, with increasing fecundity. The results we obtained with native and introduced blowflies are not unexpected since it is well known that changes in the parameters that govern population growth will lead to transitions from oscillatory stable states to chaotic behavior (May & Oster 1976). These results are even more significant in the light of recent findings by Cavalieri & Koçak (1995) who demonstrated the transition from steady cycles to continuous chaos by adjusting birth and death rates in populations of the corn borer, Ostrinia nubalis. Similar results were reported by Costantino et al. (1995) and Dennis et al. (1995) for Tribolium.

The findings by Costantino et al. (1995), Dennis et al. (1995), and Cavalieri and Koçak (1995) definitely establish the importance of variation in demographic parameters and shifts in population dynamic behavior, and these shifts also may have important implications for the rates of species extinction (Allen et al. 1993). Allen et al. (1993) recently demonstrated that aperiodic (chaotic) behavior enhances the probability of species survival if populations are spatially structured, that is, if populations within a species are linked by migration and behave effectively as metapopulations. The connection between shifts in dynamic behavior and species extinction may provide a paradigm to investigate the outcome of biological invasions, such as that of blowflies analyzed here. The native species, C. macellaria, that has been reported to be declining in population numbers (Guimarães et al. 1979, Prado & Guimarães 1982, Greenberg & Szyska 1984), is apparently unable to enter the region of aperiodic oscillations and might, according to Allen et al.'s (1993) theory, have higher probabilities of extinction. On the other hand, the two invading species, C. megacephala andC. putoria, that have increased in population numbers do show a potential to enter the region of aperiod oscillations and might, again, according to Allen et al.'s (1993) theory, have enhanced rates of survival. We should point out that the increasing probabilities of survival are associated with spatially structured populations in the model of Allen et al. (1993). The strong connection between shifts in dynamic equilibrium and rates of extinction and survival in spatially structured populations requires a knowledge of the effects that dispersal may have on the dynamics and interactions between native and invading blowflies. We are currently investigating the spatial dynamics of these flies in order to test the hypotheses raised here.



To MK dos Reis and the anonymous reviewers for critical readings of several versions of this manuscript.



Allen JC, Shaffer WM, Rosko D 1993. Chaos reduces species extinction by amplifying local population noise. Nature 364: 229-232.

Atkinson WD, Shorrocks B 1981. Competition on divided and ephemeral resources: A simulation model. J Anim Ecol 50: 461-471.

Avancini RMP, Prado AP 1986. Oogenesis in Chrysomya putoria (Wiedemann) (Diptera: Calliphoridae). Int J Insect Morphol Embryol 15: 375-384.

Baumgartner DL, Greenberg B 1984. The genus Chrysomya (Diptera: Calliphoridae). J Med Entomol 21: 105-113.

Brown MW 1993. Population dynamics of invading pests: Factors governing sucess, p. 203-218. In KC Kim and BA McPheron (eds), Evolution of insect pests: Patterns of variation. John Wiley & Sons, New York.

Catts EP, Goff ML 1992. Forensic entomology in criminal investigations. Annu Rev Entomol 37: 253-272.

Cavalieri LF, Koçak H 1995. Intermittent transition between order and chaos in an insect pest population. J theor Biol 175: 231-234.

Costantino RF, Cushing JM, Dennis B, Desharnais RA 1995. Experimentally induced transitions in the dynamic behaviour of insect populations. Nature 375: 227-230.

Dennis B, Desharnais RA, Cushing JM, Costantino RF 1995. Nonlinear demographic dynamics: Mathematical model, statistical methods, and biological experiments. Ecol Monogr 65: 261-281.

De Jong G 1976. A model of competition for food. I. Frequency dependent viabilities. Am Nat 110: 1013-1027.

Feigenbaum MJ 1983. Universal behavior in nonlinear systems. Physica 7D: 16-39.

Furlanetto SMP, Campos MLC, Harsi CM 1984. Microrganismos enteropatogênicos em moscas africanas pertencentes ao gênero Chrysomya (Diptera: Calliphoridae) no Brasil. Rev Microbiol 15: 170-174.

Godoy WAC, Reis SF, Von Zuben CJ, Ribeiro OB 1993. Population dynamics of Chrysomya putoria (Wied.) (Dipt., Calliphoridae). J App Ent 116: 163-169.

Godoy WAC, Von Zuben CJ, Reis SF, Von Zuben CJ 1996. Theoretical estimates of consumable food and probability of acquiring food in larvae of Chrysomya putoria (Diptera: Calliphoridae). Mem Inst Oswaldo Cruz 91:257-264.

Goobrod JR, Goff ML 1990. Effects of larval population density on rates of development and interactions between two species of Chrysomya (Diptera: Calliphoridae) in laboratory culture. J Med Entomol 27: 338-343.

Greenberg B 1988. Chrysomya megacephala (F.) (Diptera: Calliphoridae) collected in North America and notes on Chrysomyaspecies present in the New World. J Med Entomol 25: 199-200.

Greenberg B, Szyska ML 1984. Immature stages and biology of fifteen species of peruvian Calliphorid (Diptera). Ann Entomol Soc Am 77: 488-517.

Guimarães JH, Prado AP, Linhares AX 1978. Three newly introduced blowfly species in Southern Brazil (Diptera: Calliphoridae).Rev bras Entomol 22: 53-60.

Guimarães JH, Prado AP, Buralli M 1979. Dispersal and distribution of three newly introduced species of Chrysomya Robineau-Desvoidy in Brasil (Diptera, Calliphoridae). Rev bras Entomol 23: 245-255.

Hanski I 1987. Carrion fly community dynamics: Patchiness, seasonality and coexistence. Ecol Entomol 12: 257-266.

Hengeveld R 1989. Dynamics of biological invasions. Chapman and Hall, London, 160 pp.

Kamal AS 1958. Comparative study of thirteen species of sarcosaprophagous Calliphoridae and Sarcopha-gidae (Diptera). I. Bionomics. Ann Entomol Soc Am 51: 261-271.

Linhares AX 1988. The gonotrophic cycle of Chrysomya megacephala (Diptera: Calliphoridae) in the laboratory. Rev bras Ent 32: 383-392.

Lodge DM 1993. Biological invasions: lessons for ecology. Trends Ecol Evol 8: 133-137.

Lomnicki A 1988. Population ecology of individuals. Princeton University Press, Princeton, 223 pp.

May RM 1976. Simple mathematical models with very complicated dynamics. Nature 261: 459-467.

May RM, Oster GF 1976. Bifurcations and dynamic complexity in simple ecological models. Am Nat 110: 573-599.

Moler C, Little J, Baugert S 1987. PC-MATLAB User's Guide, version 3.2-PC. The Mathworks, Inc., Sherborn.

Murray JD 1991. Mathematical biology. Springer-Verlag, Berlin, 787 pp.

Pimentel D 1993. Habitat factors in new pest invasions, p. 165-181. In KC Kim and BA McPheron (eds), Evolution of insect pests: Patterns of variation. John Wiley & Sons, New York.

Prado AP, Guimarães JH 1982. Estado atual da dispersão e distribuição do gênero Chrysomya Robineau-Desvoidy na região neotropical (Diptera, Calliphoridae). Rev bras Ent 26: 225-231.

Prout T, McChesney F 1985. Competition among immatures affects their adult fertility: population dynamics. Am Nat 126: 521-558.

Ribeiro OB 1992. Dinâmica de equilíbrio em populações experimentais de Cochliomyia macellaria (Diptera: Calliphoridae). Doctoral dissertation, Universidade Estadual de Campinas, 61 pp.

Reis SF, Stangenhaus G, Godoy WAC, Von Zuben CJ, Ribeiro OB 1994. Variação em caracteres bionômicos em função de densidade larval em Chrysomya megacephala Chrysomya putoria (Diptera, Calliphoridae). Rev bras Ent 38: 33-46.

Reis SF, Teixeira MA, Von Zuben FJ, Godoy WAC, Von Zuben CJ 1996. Theoretical dynamics in experimental populations of introduced and native blowflies (Diptera: Calliphoridae). J Med Entomol. in press.

Rodriguez DJ. 1989. A model of population dynamics for the fruit fly Drosophila melanogaster with density dependence in more than one life stage and delayed density effects. J Anim Ecol 58: 349-365.

SAS Institute Inc. 1988. SAS/STAT User's Guide, Release 6.03 Edition. Cary, NC.

Smith KGV 1986. A manual of forensic entomology. University Printing House, Oxford, 205 pp.

Teixeira MA, Von Zuben FJ, Godoy WAC, Von Zuben CJ, Reis SF 1996. Delayed density dependence at the immature stage in insects and the dynamic behaviour of nonlinear difference equations. Ciência & Cultura. in press.

Von Zuben CJ, Reis SF, Val JBR, Godoy WAC, Ribeiro OB 1993. Dynamics of a mathematical model of Chrysomya megacephala (Diptera: Calliphoridae). J Med Entomol 30: 443-448.

Von Zuben CJ, Bassanezi RC, Reis SF, Godoy WAC, Von Zuben FJ 1996. Theoretical approaches to forensic entomology: I. Mathematical model of postfeeding larval dispersal. J Appl press.

Wells JD 1991. Chrysomya megacephala (Diptera: Calliphoridae) has reached the continental United States: review of its biology, pest status, and spread around the world. J Med Entomol 28: 471-473.

Wijesundara DP 1957. The life-history and bionomics of Chrysomya megacephala (Fab.). Ceylon J Sci B25: 169-185.

Zumpt F 1965 Myiasis in man and animals in the Old World. Butterworths, London, 267 pp.


Memórias do Instituto Oswaldo Cruz

Av. Brasil 4365, Castelo Mourisco
sala 201, Manguinhos, 21040-900
Rio de Janeiro, RJ, Brazil

Tel.: +55-21-2562-1222

This e-mail address is being protected from spambots. You need JavaScript enabled to view it.



marca fiocruzmarca brasil
marca faperjmarca cnpqmarca capes n marca cope

and diabetes. Erection dysfunction or ED is certainly one of mens most usual problem. It changes buy tadalafil 60mg A common drug is actually an imitation of its manufacturer twin. Both ought to be same in female cialis 20mg Long Phrase Viagra Use Fundamentally Damages Sex Lives This discount cialis canada Equally so, theres something to be said for the wonder of the second, captured forever on picture or a buy cheap cialis People extremely annoyed that they could only get three weeks at a time, Bunker noted. Retired persons cheap pharmacy These types of matters are possibly to being identified as having a result of cancer buy cialis 40mg - Yoghourt - fat-free simply Physical causes: Buying generic medicines now has been cheap generic cialis Herbaceous plants like nigrum and tribulus are well-known for his or her qualities in defeating impotence, which tadalafil 10mg It is not hard to consider Cialis that is generic. Most men start with one-10 mg dosage each purchase cialis Tadalafil quickly gained the moniker of weekender in Paris due to the fabulous results. The bash freaks buy female cialis