Complex Behaviour in a Discrete Coupled Logistic Model for the Symbiotic Interaction of Two Species

A symmetrical cubic discrete coupled logistic equation is proposed to model the symbiotic interaction of two isolated species. The coupling depends on the population size of both species and on a positive constant $\lambda$, named the mutual benefit. Different dynamical regimes are obtained when the mutual benefit is modified. For small $\lambda$, the species become extinct. For increasing $\lambda$, the system stabilizes in a synchronized state or oscillates in a 2 periodic orbit. For the greatest permitted values of $\lambda$, the dynamics evolves into a quasiperiodic, into a chaotic scenario or into extinction. The basins for these regimes are visualized as coloured figures on the plane. These patterns suffer different change as consequence of basins' bifurcations. The use of the critical curves let us to determine the influence of the zones with different number of first rank preimages in those bifurcation mechanisms.


DYNAMICS OF ISOLATED SPECIES: THE LOGISTIC MODEL
Imagine an island without any contact with the exterior. Living species there have no possibility to migrate looking for a new land with affordable resources. Thus, for instance, if the island has initially as inhabitants a couple of rabbits, they will reproduce exponentially. This expansion regime of the rabbits will finish for colonizing the whole island in a few generations. Hence, the island will become overpopulated. A new dynamical regime will be present now with a natural population control mechanism due to the overcrowding.
If n x represents the population after n generations, let us suppose this variable bounded in the The activation or expanding phase is controlled by the term n x µ proportional to the current population n x and to the constant growth rate µ . Resource limitations bring the system to an inhibition or contracting phase directly related with overpopulation. The term ) 1 ( can denote how far the system is of overcrowding. Therefore, if we take the product of both terms as the most simple approach to the population dynamics, the model gives account of its evolution. This is the so-called logistic map, where 4 0 < < µ in order to assure 1 0 < < n x . The continuous version of this model was originally introduced by Verhulst [Verhulst, 1845] in the nineteenth century and it has been a subject of study in the last century as a tool to be applied to the most diverse phenomenology [May, 1976] or as an object interesting to analyze by itself from a mathematical point of view [Collet & Eckmann, 1980 ;Mira, 1987].
The dynamical behavior of the logistic equation when the growth rate is modified is as follows: (i) 1 0 < < µ : The growth rate is not big enough to stabilize the population. It will drop and the specie will become extinct. (ii) 3 1 < < µ : A drastic change is obtained when µ is greater than 1. The non vanishing equilibrium between the two competing forces, reproduction on one hand and resource limitation on the other, is now possible. The population reaches, independently of its initial conditions, a fixed value that is maintained in time. (iii) 57 . 3 3 < < µ : A cascade of sudden changes provokes that the population oscillates in cycles of period n 2 , where n increases from 1, when µ is close to 3, to infinity when µ is approaching the critical value 3.57. This is named the period-doubling cascade. (iv) 82 : When the parameter moves, the system alternates between periodical behaviours with high periods on parameter interval windows and chaotic regimes for parameter values not located in intervals. The population can be not predictable although the system is deterministic. The chaotic regimes are observed for a given value of µ on sub-intervals of [0,1]. (v) 85 . 3 82 .
3 < < µ : The orbit of period 3 appears for 82 . 3 = µ after a regime where unpredictable bursts, named intermittences, have become rarer until their disappearance in the three-periodic time signal. The existence of the period 3 orbit means, such as the Sarkovskii theorem tells us, that all periods are possible for the population dynamics, although, in this case, they are not observable due to their instability. What it is observed in this range is the period-doubling cascade n 2 * 3 .
(vi) 4 85 . 3 < < µ : Chaotic behaviour with periodic windows is observed in this interval. (vii) 4 = µ : The chaotic regime is obtained on the whole interval [0,1]. This specific regime produces dynamics, which looks like random. The dynamics has lost near-all its determinism and the population evolves as a random number generator.
Therefore, there are essentially three remarkable dynamical behaviours in this system: the period doubling route to chaos around the value 57 . 3 ≈ µ [Feigenbaum, 1978], the time signal complexification by intermittence in the neighbourhood of 82 . 3 ≈ µ [Pomeau & Manneville, 1980] and the random-like dynamics for 4 = µ .

DYNAMICS OF TWO ISOLATED SYMBIOTIC SPECIES: A COUPLED LOGISTIC MODEL
Let us suppose now, under a similar scheme of expansion/contraction, that two symbiotic species ) , ( n n y x are living on the island. Each one of them evolves following a logistic-type dynamics, In this model, the symbiotic interaction between both species provokes that the growth rate ) (z µ is varying with time. In fact, it depends of the population size of the others and of a positive constant λ that we call the mutual benefit. As it is seen in the equations, we are thinking in a symmetrical interaction. Concretely, the particular dynamics of each one of the species is a logistic map whose parameter µ n is not fixed, , where λ is a real and adjustable parameter. In the following we shall write T instead of T λ as the dependence on the parameter λ is understood. Let us observe that when 0 = n y or 0 = n x , the logistic dynamics for one isolated specie is recovered. In this case the parameter λ takes the role of the parameter µ . At this point we must comment that the different choices of µ n give a wide variety of dynamical behaviours. For instance, the application of this idea produces the on-off intermittence phenomenon when µ n is chosen random [Platt et al., 1993] or the adaptation to the edge of chaos when µ n is a constant with a small time perturbation [Melby et al., 2000]. Other systems built under this mechanism are models (a), (b) and (c) presented in [López-Ruiz, 1991 ;López-Ruiz & Pérez-Garcia, 1991]. Equations (1) correspond to model (a) of that work. Model (b) has been studied in detail in , and a similar investigation of model (c) is presented in ].
In Sections 3 and 4, we study more accurately the model (1) from a dynamical point of view. In order to summarize, we explain first the dynamical behaviour of the coupled logistic system (1).
When λ is modified, it is as follows: The mutual benefit is not big enough to allow a stable coexistence of both species and they will disappear. : A sudden change is obtained when λ is greater than 0.75. Both populations are synchronized to a stable non-vanishing fixed quantity when the initial populations overcome certain critical values. If the initial species are under these limits both will become extinct (iii) The system is now bi-stable. Each one of the species oscillates out-ofphase between the same two fixed values. This is a lag-synchronized state or in other words a stable 2-period orbit. There is still in this range the possibility of extinction when the initial populations are very small or close to the overcrowding. (iv) 03 . 1 95 . 0 < < λ : The system is not anymore on a periodic orbit. It acquires a new frequency and the dynamics is now quasiperiodic. Both populations oscillate among slightly infinitely many different states. Synchronization is lost. There are in this regime periodic windows where the system becomes newly lag-synchronized. Also, for initial populations nearly zero and for 1 < λ , the species can not survive.
(v) 0843 . 1 03 . 1 < < λ : The system is now in a chaotic regime. It is characterized by a noisylike small oscillation around a synchronized state with non-periodic unpredictable bursts. Periodic oscillations can be also obtained for some particular values of the mutual benefit. Some other initial conditions are not meaningful or interpretable in this scheme because the system is going outward the square [ ] [ ] 1 , 0 1 , 0 × and evolves toward infinity. The system 'crashes'. This sudden 'damage' is interpreted as some kind of catastrophe provoking the extinction of species.
Let us remark that although the equations are formed by logistic-type components, the logistic effects have been lost and a completely new scenario emerges when they are coupled. In this case, the symbiotic interaction makes the species to reach different stable states. Depending on the mutual benefit, the system can reach the extinction, a fixed synchronized state, a bi-stable lag-synchronized configuration, an oscillating dynamics among infinitely many possible states or a chaotic regime. We must highlight in this model the phenomenon of synchronization in the periodic regime [Boccaletti et al., 2002], the transition to chaos by the Ruelle-Takens route [Eckmann, 1981] and the bursting events around a noisy-like synchronized state in the chaotic regime [Heagy et al., 1995]. All these behaviours are caused by the symbiotic coupling of the species and are not predictable from the properties of the individual logistic evolution of each one of them. Moreover, this interaction implies a mutual profit for both species. In fact, when 1 < µ one of the isolated species is extinct, but it can survive for 1 < λ if a small quantity of the other species is aggregated to the island. Hence, the symbiosis appears to be well held in this cubic model.

STABLE ATTRACTORS: SYMMETRY AND BIFURCATIONS
For the sake of clarity, firstly, we summarise the dynamical behaviour of model (1) when the mutual benefit λ is inside the interval The different parameter regions where the mapping T has stable attractors are given in the next table. The meaning of all these attractors is explained in the next sections.

Symmetry
This model has reflection symmetry P through the diagonal Note that the diagonal is T -invariant, ( ) In general, if Ω is an invariant set of T ,

( )
T Ω Ω = , so is ( ) P Ω due to the commutation property: In fact, if some bifurcation happens in the half plane below the diagonal so it is in the above half plane, and vice versa. The dynamical properties of the two halves of phase space separated by the diagonal are interconnected by the symmetry.
Also if the set Γ verifies ( ) Then the T -iteration of a reflection symmetrical set continues to keep the reflection symmetry through the diagonal. It is worth to note that the square is invariant for µ<1, but not anymore for µ>1.

Fixed Points, 2-Cycles and Closed Invariant Curves
We focus our attention on bifurcations playing an important role in the dynamics, those happening in the interval 0843 . 1 0 < < λ . In this range, there exist stable attractors for each value of λ and it will have sense to study their basins of attraction, that is, the initial populations leading to the each one of the existing final asymptotic configurations.
The restriction of T to the diagonal is a one-dimensional cubic map, which is given by the equation The restriction of the map T to the axes reduces to the logistic map ) ( , p 0 is an attractive node. For all the rest of parameter values, p 0 is a repelling node. The points ( ) , exist for every parameter value and they are unstable for every value of λ .
. < < λ , p 3 is a saddle point and p 4 is an attractive node. In this parameter interval, the whole diagonal segment between p 3 and p 4 is locus of points belonging to heteroclinic trajectories connecting the two fixed points. The point p 4 suffers a flip bifurcation when 866 . 0 2 3 ≅ = λ . It generates a stable period-two orbit p 5 6 , outside the diagonal. These points are obtained by solving the quadratic equation x . The solutions are: . , these period-2 symmetric points lose stability via a Neimark-Hopf bifurcation. The set of points p 5 6 , gives rise to a period 2 set of two stable closed invariant curves. The size of these symmetric invariant curves grow when λ increases into the interval 0 957 1 .
< < λ , and, for some values of λ , frequency locking windows are obtained.
The period two cycles on the axes appear by a period doubling bifurcation, and are found by solving the cubic equation: They have existence for λ > 3 . The solutions are: Observe that the restriction of the map T to the axes is the logistic map, so that its dynamics gives rise to the well known cyclic logistic behaviour on the axes, as explained in section 1.

Transition to Chaos
The two closed invariant curves approach the stable invariant set of the hyperbolic point p 4 on the diagonal for λ slightly larger than1 (Fig. 1a). At first sight, for λ ≈ 1 029 . , the system still seems quasi-periodic but a finer analysis reveals the fingerprints of chaotic behaviour. Effectively, a folding process takes place in the two invariant sets, cf. [Frouzakis et al., 1997], which gives rise to the phenomenon of weakly chaotic rings when the invariant set intersects itself (Fig. 1b) (cf. [Mira & al., 1996] , the tangential contact of the two symmetric invariant sets with the stable set of the saddle p 4 on the diagonal leads to the disappearance of those two weakly chaotic rings. Just after the contact, infinitely many repulsive cycles appear due to the creation of homoclinic points and a single and symmetric chaotic attractor appears (Fig. 2a). For 1 032 1 0843 .
. < < λ , this chaotic invariant set folds strongly around p 4 , and the dynamics becomes very complex (Fig. 2b). When the limit value λ = 1 084322 . is reached, the chaotic area becomes tangent to its basin boundary, the mapping iterates can escape to infinite and the attractor disappears by a contact bifurcation ( [Mira & al., 1987], chap. 5). The time behaviour of the system can be seen in Fig. 3(ac).

BASIN FRACTALIZATION
Let us now see how the different initial populations evolve toward its asymptotic stable state. This is exactly the problem of considering the basins of the different attractors of model (1). For the sake of coherence, we consider the square as the source of initial conditions having sense in our biological model, i.e., in the map T . Let us say at this point that basins constitute an interesting object of study by themselves. If a colour is given to the basin of each attractor, we obtain a coloured figure, which is a phase-plane visual representation of the asymptotic behaviour of the points of interest. The strong dependence on the parameters of this coloured figure generates a rich variety of complex patterns on the plane and gives rise to different types of basin fractalization. See, for instance, the work done by Gardini et al. [1994] and also by  in this direction. It is now our objective to analyze the parameter dependence of basin fractalization of model (1) by using the technique of critical curves.

Definitions and General Properties of Basins and Critical Curves
The set D of initial conditions that converge towards an attractor at finite distance when the number of iterations of T tends toward infinity is the basin of the attracting set at finite distance. When only one attractor exists at finite distance, D is the basin of this attractor. When several attractors at finite distance exist, D is the union of the basins of each attractor. The set D is invariant under backward connected or non-connected. A connected basin may be simply connected or multiply connected, which means connected with holes. A non-connected basin consists of a finite or infinite number of connected components, which may be simply or multiply connected. The closure of D includes also the points of the boundary D ∂ , whose sequences of images are also bounded and lay on the boundary itself. If we consider the points at infinite distance as an attractor, its basin ∞ D is the complement of the closure of D . When D is multiply connected, ∞ D is non-connected, the holes (called lakes) of D being the non-connected parts (islands) of ∞ D . Inversely, non-connected parts (islands) of D are holes of ∞ D [Mira et al., 1996].
In Sec. 3, we explained that the map (1) may possess one or two attractors at a finite distance. The points at infinity constitute the third attractor of T . Thus, if a different colour for each different basin is given we obtain a coloured pattern in the square with a maximum of two colours. In the present case, the phenomena of finite basins disappearance have their origin in the competition between the attractor at infinity (whose basin is ∞ D ) and the attractors at finite distance (whose basin is D ). When a bifurcation of D takes place, some important changes appear in the coloured figure, and, although the dynamical causes cannot be clear, the coloured pattern becomes an important visual tool to analyze the behaviour of basins.
Critical curves are an important tool used to study basin bifurcations. They were introduced by Mira in 1964 (see [Mira et al., 1996] for further details). The map T is said to be noninvertible if there exist points in state space that do not have a unique rank-one preimage under the map. Thus the state space is divided into regions, named i Z , in which points have i rank-one preimages under T. These regions are separated by the so called critical curves LC , which are the images of the curves LC −1 . If the map T is continuous and differentiable, LC −1 is the locus of points where the determinant of the Jacobean matrix of T vanishes. When initial conditions are chosen to both sides of LC , the rankone preimages appear or disappear in pairs. (See also the glossary for different technical terms used along this work).

Critical Curves and Z i -Regions of T
In our case, the map T defined in (1) is noninvertible. It has a non-unique inverse. As we know, x y x y xy x y xy x y Hence, LC −1 is independent of λ parameter and is quadratic in x and y . It can be seen that LC −1 is a curve of four branches, with two horizontal and two vertical asymptotes. The branches x + − , that multiplies the term y 2 in Eq. (2). It follows that the critical curve of rank-1, , consists of four branches. The shape of LC and LC −1 is shown in Figs. 4(a-b). LC depends on λ and separates the plane into three regions that are locus of points having 1, 3 or 5 distinct preimages of rank-1. They are named by Z i , 5 , 3 , 1 = i , respectively ( Figure  4b). Observe that the set of points with three preimages of rank-1, 3 Z , is not connected and is formed by five disconnected zones in the plane. Let us remark that LC −1 has the reflection symmetry We see in Fig. 4(b) that the four branched LC -curve divides the diagonal ∆ in five intervals. If we know the number i of preimages of rank-1 of each segment on the diagonal, the number of preimages of rank-1 of each Z i -zone of the plane is also determined. This calculation has been performed in ]. The number of rank-1 preimages of a point ( ) x x ' , ' on the diagonal can be summarised in the following table: The coordinates of the points marking the frontier between the different Z i -zones on the diagonal are: . For example, the origin p 0 is always in the Z 5 -zone. It is located into the interval limited by out of the diagonal. According to the nomenclature established in Mira et al. [1996], the map (1)

Types of Basins in T
Depending on λ , three different types of patterns are obtained in the square We proceed to present them and to explain the role played by critical curves in the bifurcations giving rise to the third basin type.

One-Colour Pattern: Extinction of Species,
In this regime, any given initial population evolves toward the extinction. The mutual benefit is small to allow the surviving of the species. Then, all initial conditions tend to zero under iteration of T . A pattern of only one colour is obtained (Fig. 5a).
If we regard the behaviour of T in the whole plane 2 R , D undergoes an interesting bifurcation consisting in the transition from a connected to a non-connected basin (Fig. 5b). It takes place when p and infinitely many small regions without connection (islands). This disaggregation is the result of infinitely many contact bifurcations, which are explained in ]. Such phenomenon can be also found in some quadratic 2 0 Z Z − maps [Mira & Rauzy, 1995]. continue to shrink to the origin, then to be extinct. The basin is a two-colour pattern (Fig. 6a). When 957 . 0 866 . 0 < < λ , 4 p bifurcates to a 2-periodic orbit and the system becomes now a lag-synchronized oscillation. The colour corresponding to this last state has gained space on the zero state in the two-colour pattern (Fig. 6b). When 1 957 . 0 < < λ , synchronization is finally lost and the system becomes quasiperiodic. Only the corners of the square

Two-Colour Pattern: Extinction or non-trivial Evolution of Species
lead to extinction in the two-colour pattern (Fig. 6c). If we regard the total basin in 2 R , D seems to be formed by the square (preimage of rank-1 of the point 0 p ), and the three triangle-shaped regions that are preimages of rank-1 of the latter component (Fig. 6d).

Two-Colour Pattern: non-trivial Evolution or Catastrophe of Species,
The extinction by exhaustion of the resources is not possible now. But a new phenomenon takes place in this range of the parameter. Some initial conditions can give rise to an evolution that surpass the boundaries of the square [ ] [ ] 1 , 0 1 , 0 × and tend to infinity. We interpret this behaviour as some kind of internal catastrophe (war, epidemics, etc.) leading also to extinction. Although we are aware of its disconcerting meaning, this would imply that an internal catastrophe can follow in this model as a consequence of the population start from some particular initial conditions. All the rest of initial conditions bring the system to a quasiperiodic state when 03 . 1 1 < < λ or a chaotic dynamical regime when 0843 . 1 03 . 1 < < λ ( Fig. 7(a-f)). Therefore, a two-colour basin is also obtained in this range of λ parameter.
In a more detailed form, the basin bifurcations happen as follows. Points ( )  p and their preimages. Due to the fact that the preimages have a finite number of accumulation points, the structure is not fractal. A similar phenomenology has been found and studied in 2 0 Z Z − maps [Mira et al., 1994]. When λ increases ( 0835 . 1 ≈ λ ), the chaotic attractor, which is limited by arcs of n LC curves, is destroyed by a contact bifurcation with its basin boundary (Fig. 7e). A new dynamical state arises. The infinite number of unstable cycles and their rank-n images belonging to the existing chaotic area before the bifurcation define a strange repulsor which manifests itself by chaotic transients (Fig. 7f). For 085 . 1 ≈ λ the basin pattern disappears definitively.
The models scattered in the literature on the three main types of population interaction, namely predator-prey situation, competition or mutualism among species, are usually stated as quadratic equations [Murray, 2002, and references therein]. In this work, we have reinterpreted a cubic twodimensional coupled logistic equation [López-Ruiz, 1991] as a discrete model to explain the evolution of two symbiotically interacting species. The symbiotic interaction between both species is population-size dependent and is controlled by a positive constant λ that we call the mutual benefit. Depending on λ , the system can reach extinction due to the small mutual benefit or the lack of resources, it can stabilize in a synchronized state or oscillates in a 2 periodic orbit for intermediate λ or it can evolve in a quasiperiodic or chaotic regime for the greatest λ . In this last scenario, there are also initial conditions leading the system to extinction. This kind of extinction has been interpreted as an internal catastrophe caused, for instance, by political decisions or by a deficient health prevision system and not by the resources to become exhausted. Another remarkable property of the model is the fact that when 1 < µ one of the isolated species is extinct, but it can survive for 1 < λ when it interacts symbiotically with one of the other species. Then, symbiosis seems to be well held in this model. Different complex colour patterns on the plane have been obtained when the mutual benefit is modified. If 75 . 0 0 < < λ , all the dynamics is attracted by the origin and a one-colour pattern is found. When 1 75 . 0 < < λ the dynamics can settle down in two possible attractors and the basins are now characterized by two colours. Finally, if 0843 . 1 1 < < λ , the two-colour basins are consequence of the two possible asymptotic states: a quasiperiodic or chaotic finite distance attractor and an additional one located at infinity.
Critical curves have been used in order to understand the basin bifurcations found in this system. Hence, common features with those present in the simplest and well-studied case of 2 0 Z Z − maps are now more evident for this map of Z Z Z Z Z  ]. Let us remark that the rich dynamics and the complex patterns produced on the plane in this model are controlled by a single parameter, in this case, the mutual benefit between the interacting species.

GLOSSARY
INVARIANT: A subset of the plane is invariant under the iteration of a map if this subset is mapped exactly onto itself.
ATTRACTING: An invariant subset of the plane is attracting if it has a neighbourhood every point of which tends asymptotically to that subset or arrives there in a finite number of iterations.
CHAOTIC AREA: An invariant subset that exhibits chaotic dynamics. A typical trajectory fills this area densely.
CHAOTIC ATTRACTOR: A chaotic area, which is attracting.
BASIN: The basin of attraction of an attracting set is the set of all points, which converge towards the attracting set.
IMMEDIATE BASIN: The largest connected part of a basin containing the attracting set.
ISLAND: Non-connected region of a basin, which does not contain the attracting set.
LAKE: Hole of a multiply connected basin. Such a hole can be an island of the basin of another attracting set.
HEADLAND: Connected component of a basin bounded by a segment of a critical curve and a segment of the immediate basin boundary of another attracting set, the preimages of which are islands.
BAY: Region bounded by a segment of a critical curve and a segment of the basin boundary, the successive images of which generate holes in this basin, which becomes multiply connected.
CONTACT BIFURCATION: Bifurcation involving the contact between the boundaries of different regions. For instance, the contact between the boundary of a chaotic attractor and the boundary of its basin of attraction or the contact between a basin boundary and a critical curve LC are examples of this kind of bifurcation.   Observe the different j Z -zones, 5 , 3 , 1 = j .