Ejaz Ahmed

Eike Ahrens

Martin Heise

Dr Michael Ruck
email: michael.ruck@tu-dresden.de

A Facile Route for the Synthesis of Polycationic Tellurium Cluster Compounds: Synthesis in Ionic Liquid Media

Keywords: Tellurium, Polycations, Cluster compounds, Ionic Liquids, Magnetic properties Tellurium, Polycations, Cluster compounds, Ionic Liquids, Magnetic properties

come

Introduction

The formation of polycationic clusters in Lewis acid is a well-known feature of tellurium chemistry. Te 4 2+ , Te 6 4+ , Te 6 2+ , Te 8 2+ , Te 7 2+ , and Te 10

2+ have been structurally characterized [1]. There are three different routes for the synthesis of tellurium polycationic clusters in literature: The classical route, performed by Gillespie and coworkers, is the oxidation of tellurium by AsF 5 or SbF 5 in liquid SO 2 using a temperature range from -196 °C to 25 °C [2]. High temperature synthesis was carried out by Corbett et al. [3] and Beck et al. [1]. Kloo and co-workers introduced the first room temperature synthesis using GaX 3 -benzene media (X = Br, Cl) [4]. _________________ Recently, we have reported a new room temperature route for the synthesis of bismuth polycationic clusters from ionic liquids and anticipated that this approach can be extended to other elements of the periodic table, e.g. tellurium [5].

The key factors of this new synthetic route are: i) Very high solubility of the elements (Bi, Te, Sb, Se, Ga, In) and their halogenides at room temperature allow reactions under mild conditions involving one-pot synthesis.

ii) Low temperature synthesis eliminates the risk of product decomposition.

iii) Easier reproducibility of the products. iv) Controllable parameters, such as solvent acidity, oxidizing/reducing agent, or halide acceptor, allow modifying the synthesis with the objective to crystallize different products.

The application of this approach in order to synthesize tellurium cluster compounds at room temperature resulted in the isolation of [Mo 2 Te 12 ]I 6 , Te 4 [AlCl 4 ] 2 , and Te 6 [WOCl 4 ] 2 . [6a]. He used a temperature gradient of 340 °C to 270 °C to get single-crystals of [Mo 2 Te 12 ]Br 6 within four weeks. Despite many attempts, the CVT synthesis of the analogous iodide failed [7]. In our studies, [Mo 2 Te 12 ]I 6 was obtained by reacting Mo, Te and I 2 in the Lewis acidic ionic liquid [BMIM]Cl/AlCl 3 at room temperature to give a dark solution within minutes. After overnight stirring and a crystallization period of only three days, singlecrystals could be isolated with a yield of 35 %.

Results and Discussions

Synthesis and Structure of [Mo 2 Te

The monoclinic cell (space group P2 1 /n) contains four formula units. The cationic complex [Mo 2 Te 12 ] 6+ has a pair of molybdenum atoms, which are connected by two bridging η 2 -Te 2 2groups (Figure 1). In addition, each molybdenum atom is coordinated by a terminal polycationic Te 4 2+ ligand in η 4 mode. As a result, both molybdenum atoms are in square-antiprismatic environment. The point symmetry of the [Mo 2 Te 12 ] 6+ complex in both compounds, i.e. in monoclinic [Mo 2 Te 12 ]I 6 and in triclinic [Mo 2 Te 12 ]Br 6 , is close to D 2h (mmm), although the site symmetry of the cluster is C i in the bromide and C s in the iodide. Moreover, there is only one independent complex present in the unit cell of [Mo 2 Te 12 ]I 6 whereas two independent complexes are found in [Mo 2 Te 12 ]Br 6 .

The distances between molybdenum and tellurium atoms of the Te 2 2-groups are slightly shorter (272.35(4) pm to 274.41(4) pm) than those between molybdenum and tellurium atoms of the Te 4 2+ rings (278.99(5) pm to 285.38(4) pm). The Te 4

2+ rings show significant distortion (distances: 276.12(5) pm to 286.85(4) pm, angles: 87.19(1) ° to 92.61(1) °) from D 4h to butterfly conformation with dihedral angles 5.3 ° and 5.9 ° (Table 1).

The complexes [(Mo 3+ ) 2 (Te 2 2-) 2 (Te 4 2+ ) 2 ] 6+ are surrounded by iodide ions, which connect them into layers (Figure 2). The shortest distance between an iodide ion and an atom of the Te 4 2+ rings is 302.20(4) pm. [Mo 2 Te 12 ]I 6 is paramagnetic with Weiss temperature θ = +5.8 K (Figure 3). An analysis based on the Curie-Weiss law gives a paramagnetic effective moment of 3.53 µ B per molybdenum atom, which is slightly lower than the spin-only value of 3.87 µ B for a d 3 configuration. In low fields some traces of ferromagnetism are found but these are very weak and most probably due to minor impurity phases.

The absence of intramolecular coupling of the magnetic cations in [Mo 2 Te 12 ]I 6 is astonishing, since the short distance of 297.16 (5) pm between the molybdenum(III) atoms suggests a chemical bond. In the present case the complex is a diradical with seventeen electrons per central atom. The electron counting scheme adds up the three valence electrons of a molybdenum(III) atom, six electrons of the π-system of the Te 4 2+ ring, and four electrons (two lone-pairs) of each Te 2 2-group (3 + 6 + 4•2 = 17).

In contrast, the closely related catenacompound MoTe 4 Br = [(Mo 3+ ) 2 (Te 2 2-) 2 (Te 2 Br -) 2 ] with a Mo-Mo bond length of 300.5(2) pm shows only weak, nearly temperature-independent paramagnetism with a magnetic moment of 0.4 µ B at 310 K, indicating a strong coupling and an eighteen electron situation [6a].

Synthesis and Structure of Te 4 [AlCl 4 ] 2

Te 4 [AlCl 4 ] 2 was reported by Corbett and coworkers through synproportionation at 250 °C [3]. Contrary to higher temperature method, we obtained Te 4 [AlCl 4 ] 2 by reacting Te and TeCl 4 in the Lewis acidic ionic liquid [BMIM]Cl/AlCl 3 at room temperature. Addition of a small amount of anhydrous NaCl is helpful for getting single crystals in two days (66 % yield). The orthorhombic compound contains planar Te 4 2+ polycations and tetrahedral AlCl 4 -anions. The atomic coordinates agree within three standard deviations with those published earlier [3].

Synthesis and Structure of Te 6 [WOCl 4 ] 2

Crystals of Te 6 [WOCl 4 ] 2 have previously been synthesized within two weeks using CVT in the temperature gradient of 230 °C to 210 °C [6b]. We carried out the oxidation of elemental tellurium with WOCl 4 in ionic liquids at room temperature. Addition of anhydrous NaCl to the dark coloured filtrate resulted in precipitation of black crystals in three days (yield: 39 %).

The crystal structure consists of boat shaped Te 6

2+ polycations and one-dimensional squarepyramidal WOCl 4 -anions. The atomic coordinates agree within six standard deviations with those published by Beck [6b].

Conclusion

We have developed an ionic-liquid-based room-temperature synthesis of homo-and heteropolyatomic cluster compounds. The significance of the present study is manifold. The redox reactions in ionic liquids provide an efficient, one-step route to the inorganic clusters. Salting out is quite helpful technique to get single crystals in few days. The additional advantages of this synthetic route are: i) Inclusion of the organic cation has not been observed yet. ii)

Ionic liquids render stability to the ionic intermediates [5]. iii)

Reactions are faster and there is no need to wait for weeks as in chemical vapor transport (CVT). Crystal growth requires only few days, and the yields can be rather high. iv)

Moreover, ionic liquids are regarded as economical and environmental friendly solvents. They could be alternatives to cancer causing organic solvents [8]. These features make the soft chemical procedure a versatile and general route towards the synthesis of other innovative materials.

Experimental Section

Chemicals: Starting materials were commercial 1methylimidazole (Merck, 99 %), 1-chlorobutane (Merck, 98.5 %), AlCl 3 (Fluka, anhydrous, 98%), Te (Aldrich, 99.999 %), WOCl 4 (Aldrich, 98 %), Mo (Aldrich, 99.95 %), I 2 (Grüssing, 99.5 %), and Cl 2 (Air Liquide, 99.8 %). TeCl 4 was synthesized from the elements [START_REF] Brauer | Handbuch der Präparativen Anorganischen Chemie[END_REF]. AlCl 3 was sublimated three times. I 2 was sublimated twice over BaO. Organic reagents were distilled before use. Because of the high moisture sensitivity of the anhydrous metal halides used in this work, all manipulations were performed under dry argon (99.999 %) atmosphere in standard Schlenk tubes.

Synthesis of [Mo 2 Te 12 ]I 6 :

The ionic liquid [BMIM]Cl/AlCl 3 ([BMIM] + : 1-n-butyl-3methylimidazolium) were prepared according to a literature procedure [START_REF] Wilkes | [END_REF]. The starting materials were added together in the Lewis acidic [BMIM]Cl/AlCl 3 ionic liquid (molar ratio = 1:1.3, volume = 1.5 ml) having molar ratio of Mo:Te:I = 1:6:3 (total mass = 355 mg). The reaction mixture was left for stirring overnight at room temperature. After three days, black, cube-shaped crystals were obtained. Yield: 35 %. The energy-dispersive Xray (EDX) analysis of a single crystal gave the atomic ratio Mo:Te:I = 10:55:35 (expected = 10:60:30). EDX is only a semiquantitative technique, which sometimes cannot give exact estimation for two closely related elements.

Synthesis of Te 4 [AlCl 4 ] 2 :

The stoichiometric amount of Te and TeCl 4 (total mass = 260 mg) were added to the Lewis acidic [BMIM]Cl/AlCl 3 (molar ratio = 1:1.3, volume = 1.5 ml) ionic liquid. The reaction mixture was left for stirring overnight at room temperature then filtered to separate unreacted material. Addition of a small amount of anhydrous NaCl into the dark violet filtrate led to the precipitation of violet, column-like crystals in two days. Yield: 66 %. EDX: Te:Al:Cl = 27:15:58 (expected = 29:14:57).

Synthesis of Te 6 [WOCl 4 ] 2 :

The same procedure was used as described for Te 4 [AlCl 4 ] 2 , starting from Te and WOCl 4 (total mass = 288 mg). After three days, black, irregular shaped crystals were obtained. Yield: 39 %. Because of the high moisture sensitivity of Te 6 [WOCl 4 ] 2 , it was not yet possible to carry out EDX analysis.

X-Ray Crystallography: Since all compounds are sensitive towards hydrolysis by moist air, crystals were selected in a glove box (concentration of O 2 < 0.1 ppm and H 2 O < 0.1 ppm) and filled in glass capillaries of 0.2 mm diameter, which were sealed by a hot filament in the argon atmosphere. Intensity data of a single-crystal of [Mo 2 Te 12 ]I 6 were recorded using a Bruker Kappa CCD diffractometer and graphite-monochromatized Mo-Kα radiation (λ = 71.073 pm) at 150 K. Data sets of single-crystals of Te 4 [AlCl 4 ] 2 and Te 6 [WOCl 4 ] 2 were collected on an imaging-plate diffraction system IPDS-I (Stoe) using graphitemonochromatized Mo-Kα radiation at room temperature. The raw data were corrected for background, polarization and the Lorentz factor. The microscopic description of the crystal shapes, which was later used in the numerical absorption corrections [START_REF]Data Reduction Program, Version 1.01[END_REF], was optimized using sets of reflections that are equivalent in the monoclinic Laue class for [Mo 2 Te 12 ]I 6 and Te 6 [WOCl 4 ] 2 , orthorhombic for Te 4 [AlCl 4 ] 2 [START_REF]Crystal Optimisation for Numerical Absorption Correction Program[END_REF]. The structures were solved with direct methods and refined using SHELX97 [13]. After refinement of the structure model of [Mo 2 Te 12 ]I 6 , we found a residual electron density of +5.08 e•10 -6 pm -3 which can be interpreted as a second orientation of the Te 4 2+ rings with weight of a few percents. This disorder could not be refined neither in a split atom model nor in a twin model. Slow cooling of the crystals did not lead to ordering. Graphics of the structures were developed using the program Diamond 3.2e [START_REF] Brandenburg | Diamond 3.2e, Crystal and Molecular Structure Visualization[END_REF]. The results of the structure determinations are summarized in Tables 2 and3. Further crystal structure data have been deposited at the Fachinformationszentrum Karlsruhe, D-76344 Eggenstein-Leopoldshafen, Germany (Fax: +49-7247-808-666, E-Mail: crysdata@fiz-karlsruhe.de). Inquiries should be accompanied by the depository number CSD-421881 ([Mo 2 Te 12 ]I 6 ), CSD-421883 (Te 4 [AlCl 4 ] 2 ), and CSD-421882 (Te 6 [WOCl 4 ] 2 ).

Magnetic Susceptibility Measurements:

The magnetic susceptibility of [Mo 2 Te 12 ]I 6 was measured with a SQUID magnetometer (MPMS-XL7, Quantum Design) between 1.8 K and 400 K of a 50 mg sample in four different magnetic fields between 0.10 and 70 kOe.

Supporting information:

Figure 1 .

 1 Figure 1. [Mo 2 Te 12 ] 6+ complex cation. The ellipsoids represent 95 % probability at 150(2) K.

Figure 2 .

 2 Figure 2. Section of the crystal structure of [Mo 2 Te 12 ]I 6. The ellipsoids represent 95 % probability at 150(2) K.

Figure 3 .

 3 Figure 3. Temperature dependence of reciprocal molar susceptibility for [Mo 2 Te 12 ]I 6 in a magnetic field of H = 35 kOe.

Figure 2 .

 2 Figure 2. Section of the crystal structure of [Mo2Te12]I6. The ellipsoids represent 95 % probability at 150(2) K. 265x187mm (96 x 96 DPI)

12 ]I 6

  [Mo 2 Te 12 ]I 6 is homeotypic to [Mo 2 Te 12 ]Br 6 containing a heteropolycationic [Mo 2 Te 12 ] 6+ cluster. [Mo 2 Te 12 ]Br 6 was synthesized by Beck through chemical vapor transport (CVT)

Table 1 .

 1 Structural parameters for Te 4 [AlCl 4 ] 2 and Te 6 [WOCl 4 ] 2 , as well as data of the temperature dependence of the reciprocal molar susceptibility for [Mo 2 Te 12 ]I 6 in different magnetic fields are available on internet http://dx.doi.org/10.1002/zaac. Selected interatomic distances (pm) in [Mo 2 Te 12 ]I 6 .

	atoms		distance	atoms		distance
	Te1-	Te2 286.85(4)	Te7-	Te8	276.81(4)
		Te4 276.97(4)		Mo1 272.73(4)
		Mo1 284.54(4)		Mo2 273.71(4)
	Te2-	Te3 276.93(4)	Te8-	Mo1 273.10(4)
		Mo1 280.18(4)		Mo2 273.49(4)
		I2	302.20(4)	Te9-	Te10 283.07(5)
	Te3-	Te4 280.66(4)		Te12 276.12(5)
		Mo1 284.54(4)		Mo2 285.38(4)
	Te4-Mo1 280.24(4)	Te10-	Te11 276.91(5)
	Te5-	Te6 276.71(4)		Mo2 279.97(5)
		Mo1 274.05(4)		I6	314.77(5)
		Mo2 274.41(4)	Te11-	Te12 281.18(5)
	Te6-	Mo1 273.22(4)		Mo2 284.43(5)
		Mo2 272.35(4)	Te12-Mo2 278.99(5)
				Mo1-Mo2 297.16(5)

Table 2 .

 2 Crystallographic data and details of the structure determinations for [Mo 2 Te 12 ]I 6 , Te 4 [AlCl 4 ] 2 , and Te 6 [WOCl 4 ] 2 .

	Compound	[Mo 2 Te 12 ]I 6	Te 4 [AlCl 4 ] 2	Te 6 [WOCl 4 ] 2
	Crystal system	monoclinic	orthorhombic	monoclinic
	Space group	P2 1 /n	Pbca	P2 1 /c
	Temperature (K)	150(2)	293(2)	293(2)
	a, b, c (pm), β (°)	1138.92(2), 1628.13(2), 1611.05(2), 105.88(1)	1072.1(1), 1412.8(1), 1183.2(2), 90	917.07(8), 1631.9(2), 3077.5(3), 117.008(7)
	V (10 6 pm 3 )	2873.42(7)	1792.1(4)	4103.4(7)
	Formula units per cell	4	4	8
	Calculated density (g•cm -3 )	5.74	3.14	4.69
	Measurement device	CCD Kappa (Bruker)	imaging plate diffractometer (IPDS-I, Stoe)
	Radiation		graphite-monochromated Mo-Kα radiation (λ = 71.073 pm)
	Measurement range	-17 ≤ h ≤ 17	-11 ≤ h ≤ 11	-9 ≤ h ≤ 9
			-27 ≤ k ≤ 26	-14 ≤ k ≤ 14	-18 ≤ k ≤ 18
			-27 ≤ l ≤ 26	-12 ≤ l ≤ 12	-34 ≤ l ≤ 34
	µ(Mo-Kα) (mm -1 )	19.25	7.70	20.60
	Measured reflections	72507	8267	23195
	Unique reflections F o > 4σ(F o )	13529, 9970	1106, 868	5789, 3911
	R(int), R(σ)	0.041, 0.039	0.079, 0.038	0.086, 0.066
	R 1 [F o > 4σ(F o )]	0.035	0.033	0.063
	R 1 (all F o )		0.060	0.050	0.103
	wR 2 (all F o	2 )	0.060	0.066	0.070
	No. of parameters	182	64	305
	Goodness of fit	1.27	1.20	1.27
	∆ρ(max, min) (e•10 -6 pm -3 )	+5.08, -4.31	+0.86, -0.76	+1.33, -1.23

Table 3 .

 3 Coordinates and equivalent isotropic displacement parameters for [Mo 2 Te 12 ]I 6 at 150(2) K. All atoms occupy Wyckoff positions 4e. U eq is defined as one third of the trace of the orthogonalized U ij tensor.

	Atom x	y	z	Ueq / pm 2
	Te1	0.31709(3) 0.37043(2) 0.17595(2) 131.6(6)
	Te2	0.36271(3) 0.38843 (2) 0.35915(2) 146.1(6)
	Te3	0.27428(3) 0.54712 (2) 0.32742(2) 155.6(6)
	Te4	0.20828(3) 0.52383(2) 0.14798(2) 155.9(6)
	Te5	-0.02771(3) 0.34316(2) 0.10664(2) 155.7(6)
	Te6	0.07955(3) 0.24667(2) 0.24851(2) 145.1(6)
	Te7	0.03788(3) 0.41154(2) 0.39261(2) 164.7(6)
	Te8	-0.07073(3) 0.50746(2) 0.25047(2) 160.4(6)
	Te9	-0.26539(3) 0.20622(2) 0.17358(2) 176.9(6)
	Te10 -0.19369(3) 0.23018(2) 0.35450(2) 207.7(7)
	Te11 -0.30688(3) 0.38207(2) 0.32504(2) 188.6(7)
	Te12 -0.35434(3) 0.36452(2) 0.14525(2) 188.0(6)
	Mo1 0.12620(3) 0.41140(2) 0.25141(2) 73.4(6)
	Mo2 -0.11637(3) 0.34246(2) 0.24880(2) 71.9(6)
	I1	0.42135(3) 0.20062(2) 0.28829(2) 171.5(6)
	I2	0.36168(3) 0.43821(2) 0.53966(2) 227.0(7)
	I3	0.16105(3) 0.44977(2) -0.04292(2) 235.9(7)
	I4	-0.33000(3) 0.29032(2) -0.04388(2) 191.8(6)
	I5	-0.07943(3) 0.05472(2) 0.29289(2) 172.1(6)
	I6	-0.14299(4) 0.32255(2) 0.53245(2) 282.7(8)

Acknowledgments. We thank Dr. W. Schnelle (Max Planck Institute for Chemical Physics of Solids, Dresden) for performing magnetic susceptibility measurements and Prof. Dr. J. Beck (University of Bonn) for helpful discussions. The authors also gratefully acknowledge the Higher Education Commission of Pakistan and the German Academic Exchange Service (DAAD) for financial support.


Mauro Mobilia
email: m.mobilia@leeds.ac.uk

Oscillatory dynamics in rock-paper-scissors games with mutations

Keywords: Population Dynamics, Fluctuations and Demographic Noise, Evolutionary Game Theory, Coevolution with Cyclic Dominance and Biodiversity, Cycles and Quasi-Cycles

We study the oscillatory dynamics in the generic three-species rock-paper-scissors games with mutations. In the mean-field limit, different behaviors are found: (a) for high mutation rate, there is a stable interior fixed point with coexistence of all species; (b) for low mutation rates, there is a region of the parameter space characterized by a limit cycle resulting from a Hopf bifurcation; (c) in the absence of mutation, there is a region where heteroclinic cycles yield oscillations of large amplitude (not robust against noise). After a discussion on the main properties of the mean-field dynamics, we investigate the stochastic version of the model within an individual-based formulation. Demographic fluctuations are therefore naturally accounted and their effects are studied using a diffusion theory complemented by numerical simulations. It is thus shown that persistent erratic oscillations (quasi-cycles) of large amplitude emerge from a noise-induced resonance phenomenon. We also analytically and numerically compute the average escape time necessary to reach a (quasi-)cycle on which the system oscillates at a given amplitude.

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

Understanding the mechanisms allowing co-evolution and the maintenance of biodiversity is a key issue in theoretical biology and ecology [START_REF] May | Stability and Complexity in Model Ecosystems[END_REF]; [START_REF] Michod | Darwinian Dynamics[END_REF]; [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeography[END_REF]; [START_REF] Haken | Synergetics[END_REF]; [START_REF] Murray | Mathematical Biology[END_REF]; [START_REF] Neal | Introduction to Population Biology[END_REF]). In this context, cyclic dominance has been identified as a potential mechanism that helps promote species diversity (Durrett andLevin (1994, 1997); [START_REF] Gilg | Cyclic dynamics in a simple vertebrate predator-prey-community[END_REF]; [START_REF] Kerr | Local dispersal promotes biodiversity in a real-life game of rock-paper-scissors[END_REF]; [START_REF] Czárán | Chemical warfare between microbes promotes biodiversity[END_REF])

and is naturally investigated in the framework of evolutionary game theory (EGT) [START_REF] Smith | Evolution and the Theory of Games[END_REF]; [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF]; [START_REF] Nowak | Evolutionary Dynamics[END_REF]; [START_REF] Szabó | Evolutionary games on graphs[END_REF]). The emergence of oscillatory (or quasioscillatory) behavior is one of the most appealing and debated phenomena that often characterizes co-evolution in population dynamics [START_REF] Sinervo | The rock-scissors-paper game and the evolution of alternative male strategies[END_REF]; Zamudio (2000); [START_REF] Kerr | Local dispersal promotes biodiversity in a real-life game of rock-paper-scissors[END_REF]; [START_REF] Kirkup | Antibiotic-mediated antagonism leads to a bacterial game of rock-paper-scissors in vivo[END_REF]; [START_REF] Dawkins | The Selfish Gene[END_REF]; Bastolla et al. (2002); [START_REF] Hauert | Volunteering as Red Queen mechanism for cooperation in public goods games[END_REF]; [START_REF] Mobilia | Phase Transitions and Spatio-Temporal Fluctuations in Stochastic Lattice Lotka-Volterra Models[END_REF]). Oscillatory dynamics has notably been observed in predator-prey and host-pathogen systems [START_REF] Anderson | The Invasion, Persistence and Spread of Infectious Diseases within Animal and Plant Communities[END_REF]May (1986, 1991); [START_REF] Berryman | Population Cycles[END_REF]; Turchin ( 2003)), as well as in genetic networks [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF]). Here, we study the oscillatory dynamics of rock-paper-scissors games with mutations and demonstrate that this favors the longlasting co-evolution of all species. Rock-paper-scissors games (RPS) -in which rock crushes scissors, scissors cut paper, and paper wraps rock [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF]; [START_REF] Nowak | Evolutionary Dynamics[END_REF]; [START_REF] Szabó | Evolutionary games on graphs[END_REF])have emerged as paradigmatic mathematical models in EGT to describe the cyclic competition in ecosystems (Tainaka (1994); [START_REF] May | Nonlinear aspects of competition between species[END_REF]; Szabó and Szolnoki (2002); [START_REF] Reichenbach | Coexistence versus extinction in the stochastic cyclic Lotka-Volterra model[END_REF]Reichenbach et al. ( , 2007aReichenbach et al. ( , 2008b)); [START_REF] Perc | Cyclical interactions with alliance-specific heterogeneous invasion rates[END_REF]; [START_REF] Claussen | Cyclic Dominance and Biodiversity in Well-Mixed Populations[END_REF]; [START_REF] Peltomäki | Three-and four-state rock-paper-scissors games with diffusion[END_REF]; [START_REF] Berr | Zero-one survival behavior of cyclically competing species[END_REF]). In addition to their theoretical relevance, and in spite of their apparent simplicity, it has been argued that the RPS model and the like can help understand the co-evolutionary dynamics of different biological systems, such as the cyclic dominance observed in some communities of lizards [START_REF] Sinervo | The rock-scissors-paper game and the evolution of alternative male strategies[END_REF]; Zamudio (2000)).

Another popular example is the cyclic competition between three strains of E.coli [START_REF] Kerr | Local dispersal promotes biodiversity in a real-life game of rock-paper-scissors[END_REF]). In this case, it has been found that cyclic dominance yields coexistence of all species in a spatial setting, while in a well-mixed (homogeneous) environment two species go extinct after a short transient. Naturally therefore, the mathematical properties of the RPS model and its variants have recently attracted much interest. In an EGT setting, the RPS dynamics is classically described in terms of the replicator equations (REs) that is a set of deterministic rate equations (see below) [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF]; [START_REF] Nowak | Evolutionary Dynamics[END_REF]; [START_REF] Szabó | Evolutionary games on graphs[END_REF]).

The latter essentially predict two types of behavior for the RPS dynamics: there is either the stable coexistence of all species, or oscillations of large amplitudes with each species almost taking over the entire system in turn and then suddenly almost going extinct. As population dynamics always involves a finite (yet, often large) number of discrete entities, it is necessary to take the effects of noise into account and understand its influence on the ensuing nonlinear dynamics [START_REF] Nisbet | Modelling fluctuating populations[END_REF]; [START_REF] Ewens | Mathematical Population Genetics[END_REF]). In the context of population models, demographic stochasticity (i.e. intrinsic noise) is naturally accounted by adopting an 'individual-based' modeling. The stochastic dynamics is thus implemented in terms of random birth and death events, like in (1994); Reichenbach et al. (2007aReichenbach et al. ( ,b, 2008b)); [START_REF] Peltomäki | Three-and four-state rock-paper-scissors games with diffusion[END_REF]). On the other hand, for systems with homogeneous (well-mixed) populations, it is well established that intrinsic noise drastically alters the replicator dynamics and recent studies have focused on the mean time and probability of extinction of two species [START_REF] Reichenbach | Coexistence versus extinction in the stochastic cyclic Lotka-Volterra model[END_REF]; [START_REF] Claussen | Cyclic Dominance and Biodiversity in Well-Mixed Populations[END_REF]; [START_REF] Berr | Zero-one survival behavior of cyclically competing species[END_REF]).

In addition to selection and reproduction, a third basic evolutionary mechanism is mutation. The latter is generally regarded as providing the advantageous traits that survive and multiply in offsprings [START_REF] Neal | Introduction to Population Biology[END_REF]; [START_REF] Michod | Darwinian Dynamics[END_REF]). Also, from a behavioral perspective, it is recognized that agents never behave perfectly rationally and it is sensible to assume that individuals can switch their strategies and therefore undergo mutations [START_REF] Antal | Strategy abundance in 2x2 games for arbitrary mutation rates[END_REF]; [START_REF] Szabó | Evolutionary games on graphs[END_REF]). While selection and reproduction are present in the RPS models considered in Refs. (Tainaka (1994); Szabó and Szolnoki (2002); [START_REF] Reichenbach | Coexistence versus extinction in the stochastic cyclic Lotka-Volterra model[END_REF]Reichenbach et al. ( , 2007a)); [START_REF] Perc | Cyclical interactions with alliance-specific heterogeneous invasion rates[END_REF]; [START_REF] Claussen | Cyclic Dominance and Biodiversity in Well-Mixed Populations[END_REF]; Reichenbach et al. (2007bReichenbach et al. ( , 2008b)); [START_REF] Peltomäki | Three-and four-state rock-paper-scissors games with diffusion[END_REF]; [START_REF] Berr | Zero-one survival behavior of cyclically competing species[END_REF]), the latter do not account for mutations. Here, we study the oscillatory dynamics of RPS games with mutations and show that the possibility to switch from one strategy (species) to another with a small transition rate yields novel oscillatory dynamics and favors long-term coexistence. For population of finite size, we discuss the effect of demographic noise on the dynamics and show that it induces persistent (quasi-)cyclic behavior characterized by sustained erratic oscillations of non-vanishing amplitude.

The remainder of the paper is organized along the following lines. In the next section we introduce the generic RPS model with mutations (RPSM) and describe its dynamics in the mean-field limit. The (generalized) REs are studied in Section 3, where the bifurcation diagram is obtained (Sec. 3.1) and a new region of the parameter space characterized by a Hopf bifurcation is identified. The main properties of the resulting limit cycle are briefly discussed in Sec. 3.2. Section 4 is dedicated to the stochastic dynamics of the model RPSM in terms of an individual-based formulation. In particular, we show that demographic noise can cause quasi-cyclic behavior with persistent oscillations of large amplitude (Sec. 4.1), and allow to escape -after a characteristic time that is computed -from the coexistence fixed point and reach a given (quasi-)cycle (Sec. 4.2). In the final section, we summarize and discuss our results.

Rock-Paper-Scissors with Mutations

In their essence, all variants of the RPS game aim at describing the co-evolutionary dynamics of three species, say A, B and C, in cyclic competition. In this setting, as in the children's game, 'rock' (here species A) crushes the 'scissors' (species B), and 'paper' (species C) wraps the 'rock', and 'scissors' cut the 'paper'. We therefore

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

say that A dominates over C, which outcompetes B, which outgrows A and thus closes the cycle. In EGT, the interactions are specified in terms of a payoff matrix P. Generically, the cyclic dominance of RPS games is captured by the following payoff matrix [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF]; [START_REF] Nowak | Evolutionary Dynamics[END_REF]; [START_REF] Szabó | Evolutionary games on graphs[END_REF]),

where > 0:

P = vs Rock (A) Paper (B) Scissors (C) Rock (A) 0 - 1 Paper (B) 1 0 - Scissors (C) - 1 0
According to this matrix, when a pair of A and B players interacts, the former gets a negative payoffwhile the latter gets a payoff 1. In this case, A is dominated by B and its loss is less than B's gain when 0 < < 1, whereas B's gain is higher than A's loss when > 1. In the same way, C's dominate over B's and the latter prevail against A's, and thus the model exhibits cyclic dominance. In all pairwise interactions, the dominant individual gets a payoff 1, whereas the dominated one obtain a negative payoff -. Therefore, the parameter allows to introduce an asymmetry (when = 1) in the interactions. When = 1, one of the player loses what the other gains and this perfect balance corresponds to a zero-sum game.

In EGT, one often considers homogeneous (well-mixed) populations of N individuals, with N → ∞. In this mean-field limit, the dynamics is usually specified in terms of the REs for the densities (or relative abundances) a(t), b(t) and c(t) of species A, B and C, respectively [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF]; [START_REF] Nowak | Evolutionary Dynamics[END_REF]; [START_REF] Szabó | Evolutionary games on graphs[END_REF]). Introducing the vector s(t), whose components s i , with i ∈ (A, B, C), are

s A ≡ a(t), s B ≡ b(t) and s C ≡ c(t), the REs read ṡi = s i [(Ps) i -s.Ps] = s i [π i -π],
where the dot stands for the time derivative. The important notion of average payoff (per individual) of species i, π i , has been introduced in terms of the payoff matrix as a linear function of the relative abundances: π i ≡ (Ps) i , whereas π = s.Ps = i s i π i denotes the population's mean payoff [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF]; [START_REF] Nowak | Evolutionary Dynamics[END_REF]; [START_REF] Szabó | Evolutionary games on graphs[END_REF]). In addition to the above processes of selection/reproduction, we introduce a third evolutionary mechanism that allows each individual to mutate from one species to another with rate μ:

A μ -→ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ B C , B μ -→ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ A C , C μ -→ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ A B .
(1)

In this setting, the natural generalization of the replicator equations for the model under consideration is Here, μ c = 0.0125 (see text). Left: In the absence of mutation rate (μ = 0, λ = 0.0375), the fixed point s * is unstable and the density (brown, dashed) jumps from 0 to 1. For very low values of the mutation rate, here μ = 0.001 (λ = 0.0345), s * is also unstable and a(t) oscillates regularly (magenta, solid) with an amplitude that approaches the extreme values 0 and 1, but with a much shorter period than in the case μ = 0. Right: For μ = 0.01 (λ = 0.0075), the fixed point s * is still unstable and a(t) oscillates regularly about a * = 1/3 (red, solid), with a period T ≈ 10. When μ = 0.02 > μ c (λ = -0.0225), the densities exhibit exponentially damped oscillations and converge towards the fixed point value 1/3 (blue, dashed).

ṡi = s i [π i -π] + μ(1 -3s i ), with π A = c -b, π B = a -c, π C = b -

The Rate Equations

Here, we consider the dynamics of the model RPSM in the mean-field limit, where N → ∞. When the population is well-mixed (homogeneous) and every pair of random individuals has the same probability to interact, demographic noise can be neglected and the system's dynamics is aptly described by the rate equations

ȧ = a [c -b -(1 -) {ab + bc + ac}] + μ(1 -3a) ḃ = b [a -c -(1 -) {ab + bc + ac}] + μ(1 -3b) ( 2 ) ċ = c [b -a -(1 -) {ab + bc + ac}] + μ(1 -3c)
As the system is comprised of three species, the sum of the population density is conserved, i.e. a(t)+b(t)+c(t) = 1. Clearly, this constant of motion allows, say, to set c(t) = 1a(t)b(t) and reduces (2) to a system with only two variables (here, a(t) and b(t)). Thereafter, the properties of these equations are studied and it will be showed that mutations yield a new region in the parameter space characterized by a stable limit cycle.

It is known that in the absence of mutations (μ = 0) the REs admit one interior fixed point,

s * = (a * , b * , c * ) =
(1/3, 1/3, 1/3), and three absorbing states ((1, 0, 0), (0, 1, 0) and (0, 0, 1)). The resulting cyclic dynamics gives rise to three kinds of behavior [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF]; [START_REF] Nowak | Evolutionary Dynamics[END_REF]; [START_REF] Szabó | Evolutionary games on graphs[END_REF]): (i) When < 1, s * is the system's only attractor, it is globally stable and the trajectories in the phase portrait spiral towards it (s * is a focus). (ii) When > 1, the interior rest point s * becomes unstable and the flows in the phase portrait form a heteroclinic cycle connecting each of the absorbing fixed points (saddles) at the boundary of the phase portrait. Clearly, as any fluctuations cause the rapid extinction of two species, the resulting oscillations > 1 and μ ≥ 0 for ≤ 1, s * is a stable focus. Below the line μ c ( ), in region (b) where 0 < μ < μ c , the system undergoes a Hopf bifurcation (at μ = μ c ) and the dynamics is characterized by a limit cycle and stable oscillations. In region (c), where μ = 0 (no mutations) and > 1 (dashed line), the REs (2) yield heteroclinic cycles along with periodic oscillations unstable against demographic fluctuations. are non-robust [START_REF] May | Nonlinear aspects of competition between species[END_REF]; [START_REF] Durrett | The importance of being discrete (and spatial)[END_REF]). (iii) When = 1, π vanishes (zero-sum game) and the interior rest point s * is marginally stable (center). In this case, the quantity a(t)b(t)c(t) is a constant of motion and the trajectories in the phase portrait form closed orbits around s * , set by the initial conditions and not robust against noise (see e.g. [START_REF] Reichenbach | Coexistence versus extinction in the stochastic cyclic Lotka-Volterra model[END_REF])).

In the presence of mutations, solving Eqs. ( 2) with ȧ = ḃ = ċ = 0, one finds that s * = (1/3, 1/3, 1/3) is the only (interior) fixed point of the system. The absence of absorbing points when μ > 0 indicates that the model under consideration does not yield heteroclinic cycles. In fact, as shown below, in the presence of small mutation rate the heteroclinic cycles are replaced by stable cycles resulting from a Hopf bifurcation (region (b) in Fig. 2) and yield persistent oscillatory dynamics of all species, which is a feature of biological relevance.

Linear stability analysis

To investigate the properties of REs (2), it is useful to introduce the variables x = (x A , x B ):

x A = a -a * = a - 1 3 and x B = b -b * = b - 1 3 , ( 3
)
which measure the deviations from the interior fixed point s * . We notice that

x C = c -c * = -(x A + x B ).
As a first step to gain some insight into the dynamics described by (2), it is useful to perform a linear stability analysis. To linear order in terms of x = (x A , x B )), the REs (2) can be rewritten as ẋ = A(s * )x, where A(s * ) is the Jacobian matrix at s * whose (complex conjugate) eigenvalues which are

λ ± = λ ± i ω 0 , with λ = -1 6 -3μ and ω 0 = 1 + 2 √ 3 . ( 4
)

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

General results therefore guarantee that s * is (globally) stable when λ < 0, i.e. for 18μ > -1, whereas s * is unstable when λ > 0 and a Hopf bifurcation occurs at the critical value μ = μ c ( ) and > 1 is kept fixed, with

μ c ≡ -1 18 > 0, (5)
The bifurcation diagram of the model is thus summarized in Fig. 2. Such a diagram clearly illustrates the effect of mutations and is characterized by three regions. In region (b), where 0 < μ < μ c , s * is unstable and there is a Hopf bifurcation at μ = μ c (as demonstrated below) which leads to the coexistence of all species with stable oscillations of their densities (right panel of Fig. 1, red/solid). On the other hand, in region (a),

where μ >

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ μ c , for > 1 0 , for ≤ 1
, i.e. when λ < 0, the densities are exponentially damped, with an amplitude A amp (t) ≈ e -|λ|t that vanishes near s * (right panel of Fig. 1, blue/dashed). In the absence of mutations and when > 1 (region (c) in Fig. 2), one recovers the scenario characterized by heteroclinic cycles (left panel of Fig. 1, brown/dashed). In this case, the density of each species oscillates, with a large period, between the extreme values 0 (absence of a given species) and 1 (presence of only one species). These oscillations are not robust against demographic noise: chance fluctuations will unavoidably and quickly drive the system into an absorbing state. When λ = 0, the fixed point s * is a center and its stability is determined by nonlinear terms.

Below, we shall consider the dynamics to third order about s * and show that it is stable when λ = 0 and μ > 0.

Only in the marginal case λ = μ = 0, s * is marginally stable (see e.g. [START_REF] Reichenbach | Coexistence versus extinction in the stochastic cyclic Lotka-Volterra model[END_REF])).

Hopf bifurcation and limit cycle

As qualitatively discussed above, at mean-field level, the most interesting effect of mutations is the emergence of robust oscillatory behavior (stable limit cycle) resulting from a Hopf bifurcation occurring at low mutation rates (right panel of Fig. 1, red/solid). To study the properties of the resulting limit cycle and take advantage of the system's symmetry around s * , it is convenient to perform the linear transformation x → Sx ≡ y, where 4) and ( 7).

S = ⎛ ⎜ ⎝ 1 1 2 0 √ 3 2 ⎞ ⎟ ⎠.
We notice that r ∞ increases when is raised and μ is lowered.

where

α = 18ω 0 (1 + 2 √ 3ω 0 ) 7(1 + 2 ) + (13 -9μ) + 9μ(1 + 9μ) (7) β = 1 -- 6λ(1 + 2 √ 3ω 0 ) 7(1 + 2 ) + (13 -9μ) + 9μ(1 + 9μ) .
One readily recognizes that (6) yield a supercritical Hopf bifurcation when λ > 0 (for 0 < μ < μ c ( )). In fact, the parameter β is the first Lyapunov coefficient and is negative when λ > 0 and the mutation rate μ is small (i.e. μ/ 1) (Wiggins (1990); Grimshaw (1993); [START_REF] Strogatz | Nonlinear dynamics and chaos[END_REF]).

It is insightful to solve the radial component of the Eq. ( 6), whose solution is

r(t) = r(0) e λt √ 1-r(0) β λ (e 2λt -1)
. Thus, when λ > 0, the interior fixed point s * is unstable and the long-time behavior is

r(t) r ∞ 1 -e -2λt 2 λ-r 2 (0)|β| r 2 (0)|β|
, where

r ∞ = λ |β| . ( 8
)
The dependence of r ∞ on the parameters and μ is illustrated in Figure 3, where it is shown that r ∞ increases when is raised and μ lowered. Thus, it follows from Eqs. (6,8) that the REs (2) indeed yield a stable limit

cycle σ(t) = (ā(t), b(t)) (with c(t) = 1 -ā(t) -b(t))
of radius r ∞ , period T and frequency ω, where

T = 2π ω , with ω = ω 0 - λα |β| (9)
We notice that the linear terms contribute to the frequency ω through ω 0 (natural frequency), while nonlinearity gives rise to an additional contribution (-λα/|β|). According to (6), the limit cycle σ(t) is approached As they result from a third-order expansion, the predictions of ( 6) are accurate for mutation rates close to the critical value μ c (with fixed), i.e. for "small" positive values of λ and r ∞ as illustrated in Fig. 4.

exponentially fast in time, r(t) -r ∞ ∝ e -2λt .
To conclude this section, we also consider the case where λ < 0 and notice that Eqs. ( 6) thus yield r(t) t→∞ -→ e -|λ|t . This implies that s * is approached in an oscillatory manner with an exponentially damped amplitude, namely x A (t) ∝ e -|λ|t cos (ω 0 t)sin (ω 0 t)/ √ 3 and x B (t) ∝ e -|λ|t sin (ω 0 t). Furthermore, in the marginal case λ = 0 (with μ > 0), s * is still stable but approached in a slow oscillatory manner. In fact, it follows from (6) that in this case r(t) ∼ t -1/2 , which implies that the oscillations are damped algebraically as t -1/2 . Indeed,

for λ = 0 one finds x A (t) t→∞ -→ cos (ω 0 t) -1 √ 3 sin (ω 0 t) /(6 √ μt) and x B (t) t→∞ -→ sin (ω 0 t)/(3 √ 3μt).

The Stochastic Nonlinear Dynamics

We so far have focused on the deterministic description of the model in the mean-field limit. However, as virtually all real systems are of finite size, N < ∞, and made of discrete entities, they are influenced by stochastic fluctuations. Here, adopting an individual-based approach, we discuss how demographic noise alters the deterministic properties of the model RPSM.

Before mathematically studying the model RPSM in the presence of demographic noise, it is worth gaining some insight into its stochastic dynamics. As explained below, the latter has been implemented by a birthdeath process simulated according to the Gillespie algorithm [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF][START_REF] Gillespie | Exact simulations of coupled chemical reactions[END_REF]). Some characteristic results are reported in Figs. 567. The behavior corresponding to the regime (a) of the parameter space (see Fig. 2) is illustrated in Fig. 5, where demographic noise causes erratic oscillations forming "(phase-forgetting) quasi-cycles" [START_REF] Nisbet | Modelling fluctuating populations[END_REF]) in the phase portrait (see also Fig. 7, cyan/light grey trajectory).

When the population size N is large and sample-averaged over many replicates, the amplitude of the stochastic oscillations decreases and the predictions of the rate equations are approached (thick curves in Fig. 5), but fully recovered only in the mean-field limit N → ∞. The stochastic dynamics in the region (b) of the parameter space, where s * is unstable, is illustrated in Fig. 6 and is characterized by persistent oscillations with erratic, but non-vanishing, amplitude and phase. As a result, instead of a perfectly closed orbit the trajectories in the phase portrait form a perturbed limit cycle, also called "phase-remembering quasi-cycle" [START_REF] Nisbet | Modelling fluctuating populations[END_REF]), as shown in Fig. 7 (red/dark grey trajectory). ). We will also see how noise affects the system's attractors that can be regarded as minima of an effective potential well from which it takes an enormous amount of time to escape [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]; [START_REF] Kubo | Fluctuation and relaxation of macrovariables[END_REF]; [START_REF] Dykman | Large fluctuations and optimal paths in chemical kinetics[END_REF]; Volovik et al. ( 2009)).

To describe the stochastic nonlinear dynamics of the system, we adopt an individual-based approach and consider an urn model [START_REF] Feller | An Introduction to Probability Theory and its Application[END_REF]) comprising a total of and the offspring replaces a randomly chosen individual of the other species (B in this example). The fMP therefore conserves the total size N of the population and can be regarded as a random-walk with hopping rates depending on each species average payoff (fitness). While there are various ways of implementing the dynamics of the fMP, we here consider that in the absence of mutations and large (yet finite) population size the transition T B→A from A to B, say, is given by [1 + (π Aπ)] ab. This expression comprises a contribution proportional to the average payoff differences (π Aπ), that accounts for selection, supplemented by a constant (set to 1) accounting for the background random noise [START_REF] Nowak | Emergence of cooperation and evolutionary stability in finite populations[END_REF]; [START_REF] Nowak | Evolutionary Dynamics[END_REF]). The multiplying factor ab encodes the probability (when N 1) that two individual of species A and B interact. The effect of mutations is thus taken into account by adding a linear term, which finally leads to the following transition rate: T B→A (s) = (1 + {π A -π}) ab + μb. More generally, we consider the following transition rates:

N = N A + N B + N C individuals, N A are of species A, N B of
T i→j (s) = (1 + {π j -π}) s i s j + μs i , ( 10
) with i = j ∈ (A, B, C) and π A = c -b, π B = a -c, π C = b -a and π = (1 -)(ab + bc + ac).
The stochastic description of the system is encoded in the probability P (N A , N B ; t) of having N A and N B individuals of species A and B, respectively, at time t. The quantity P (N A , N B ; t) obeys the master equation associated with the above fMP (see, e.g., [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]; [START_REF] Ewens | Mathematical Population Genetics[END_REF])) and specified by the transition rates (10). The Markovian process defined by the transition rates (10) has been simulated using the Gillespie algorithm [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF][START_REF] Gillespie | Exact simulations of coupled chemical reactions[END_REF]) which was initially introduced to simulate chemical systems through their "microscopic reactions". In the vicinity of s * and for large (yet finite) population sizes, the master equation of the above fMP can be aptly described within a generalized diffusion approximation for the probability density P (x, t), where x = (x A , x B ) and x j ≡ When μ = 0.008, the interior fixed point s * is unstable and each trajectory in the phase portrait forms a "phaseremembering quasi-cycle"of erratic non-vanishing amplitude and phase orbiting around σ. When μ = 0.05, the flow spirals towards the stable interior fixed point s * and erratically wanders in its vicinity forming a "phase-forgetting quasi-cycle" (see we expect fluctuations of order N -1/2 , see below). The temporal development of P (x A , x B ; t) is thus described by a forward Kolmogorov equation (FKE), or Fokker-Planck equation, resulting from a size-expansion of the master equation (see, e.g., [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]; [START_REF] Ewens | Mathematical Population Genetics[END_REF]; [START_REF] Nisbet | Modelling fluctuating populations[END_REF]; Traulsen et al. ( 2005))).

Performing a (linear) van Kampen expansion about the fixed point s * in the continuum limit (assuming that N is large but finite), the standard procedure (see, e.g., [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]; [START_REF] Risken | The Fokker-Planck Equation[END_REF]; [START_REF] Reichenbach | Coexistence versus extinction in the stochastic cyclic Lotka-Volterra model[END_REF])) leads to the following FKE:

∂ t P (x, t) = -∂ xi [x j A ij (s * )P (x, t)] + 1 2 B ij (s * )∂ xi ∂ xj P (x, t), ( 11
)
where we have adopted the summation convention on repeated indices i, j ∈ (A, B). In Eq. ( 11) the drift terms (A(s * ) x) i are obtained from the Jacobian matrix A of the REs (2) at s * (see Sec. 3.1), while the diffusion matrix B is defined by (see, e.g., [START_REF] Claussen | Cyclic Dominance and Biodiversity in Well-Mixed Populations[END_REF]))

B AA (s * ) = T B→A (s * ) + T A→B (s * ) + T C→A (s * ) + T A→C (s * ) = 4(1 + 3μ) 9N B BB (s * ) = T A→B (s * ) + T B→A (s * ) + T C→B (s * ) + T B→C (s * ) = 4(1 + 3μ) 9N (12) B AB (s * ) = B BA (s * ) = -T A→B (s * ) + T B→A (s * ) = - 2(1 + 3μ) 9N .
As well known, the FKE ( 11) is equivalent to the following set of linear stochastic (Langevin) differential equations with white noise vector η = (η A , η B ) [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]; [START_REF] Risken | The Fokker-Planck Equation[END_REF]):

ẋA = A AA (s * ) x A + A AB (s * ) x B + η A ẋB = A BA (s * ) x A + A BB (s * ) x B + η B , (13)
with η i = 0 and covariance matrix η i (t)η j (t ) = B ij (s * ) δ(tt ), where i, j ∈ (A, B) and ... denotes the ensemble average over a large number of replicates. As anticipated, it follows from (12,13) that the noise strength is ∝ N -1/2 and yields fluctuations of order O(N -1/2 ) around the deterministic values of x A and x B .

Below, the discussion of stochastic effects on the evolutionary dynamics of the model are centred two aspects: (i) we show how demographic noise leads to the emergence of quasi-cycles; (ii) and we compute the average time for a system with the same initial density each species to "escape" from the interior rest point and reach a given cycle.

Quasi-cycles and noise-induced resonance amplification

In this section, we are especially interested in the stochastic dynamics in region (a) of the phase parameter, where λ < 0 and s * is a stable fixed point. Here, we investigate one of the most intriguing effects of intrinsic noise, which yields persistent erratic oscillations around s * (see Fig. 5) and therefore considerably alters the predictions of the REs (2). In fact, while one could naively expect only small corrections (of order N -1/2 , see Eqs. (12,13))

to the deterministic predictions, it has been suggested that demographic noise can be sufficient to perturb the stationary state predicted by the mean-field analysis and to produce persistent erratic oscillations (see e.g. [START_REF] Bartlett | Measles periodicity and community size[END_REF][START_REF] Bartlett | The critical community size for measles in the United States[END_REF])). This is indeed illustrated by the numerical simulations reported in Fig. 5. It has recently been shown that such a behavior, often referred to as quasi-cyclic [START_REF] Nisbet | Modelling fluctuating populations[END_REF]), To analytically demonstrate this phenomenon, we shall compute the power-spectrum P (Ω) of the system. The power-spectrum is one of the most useful tools to look for oscillations in noisy data and is related to the auto-correlation functions of (stationary) systems [START_REF] Nisbet | Modelling fluctuating populations[END_REF]; [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]). Introducing

x(Ω) ≡ 1 2π ∞ -∞ dt e -iΩt
x(t), the Fourier transform of x(t), the power spectrum is given by the following ensemble average:

P (Ω) = | x(Ω)| 2 = 2 | x A (Ω)| 2 = 2 | x B (Ω)| 2 .
Taking the Fourier transform of Eqs. ( 13) and then averaging the square modulus of x(Ω) one finds:

P (Ω) = 8(1 + 3μ) 9N Ω 2 0 + Ω 2 (Ω 2 -Ω 2 0 ) 2 + (2λΩ) 2 , with 9Ω 2 0 ≡ 1 + 2 √ 3 ω 0 + 9μ(9μ + 1 -). ( 14
)
For low mutation rate μ (i.e. when μ/ 1), Ω2 0 > 2λ 2 and the power spectrum P (Ω) has a single peak (for Ω > 0) at the characteristic frequency

Ω * = Ω 0 ⎛ ⎝ 2 1 - λ Ω 0 2 -1 ⎞ ⎠ 1/2 , ( 15
)
as illustrated in Fig. 8. When the mutation rate μ is low, the value of Ω * is very close to (yet different from) ω 0 and Ω 0 (4, 14) 2 . When (λ/Ω 0 ) 2 1, the above expression simplifies to give Ω * = Ω 0 (1 + O(λ 2 /Ω 2 0 )). It is also worth noticing that at the onset of the Hopf bifurcation, i.e. when λ = 0 (with μ > 0), the denominator of P (Ω)

vanishes for Ω = Ω 0 and in this case the peak in the power spectrum is replaced by a pole at Ω = Ω * = Ω 0 .

While on general grounds (law of large numbers [START_REF] Feller | An Introduction to Probability Theory and its Application[END_REF])) one expects oscillations of the order N -1/2 , the presence of a peak in the power spectrum corresponds to a noise-induced amplification due to a resonance effect. In fact, there is a resonant behavior when there exists a particular frequency for which the denominator of the power spectrum P (Ω) is small. Here, the denominator of the power spectrum takes its minimal value at

Ω res = Ω 2 0 -2λ 2 , ( 16
)
which simplifies when (λ/Ω 0 ) 2 1 and also gives Ω res = Ω 0 (1 + O(λ 2 /Ω 2 0 )). In such a regime, where Ω * Ω res Ω 0 , there is a very suggestive explanation of the noise-induced resonant amplification by comparison with a simple mechanical system [START_REF] Mckane | Predator-prey cycles from resonant amplification of demographic stochasticity[END_REF]). Indeed, one can consider a linear damped harmonic oscillator of natural frequency Ω 0 , with damping constant ζ, driven by a force oscillating with frequency Ω. One can thus show that such a mechanical system oscillates at frequency Ω, with amplitude A amp (Ω) ∝

(Ω 2 -Ω 2 0 ) + (ζΩ) 2 -1 and yields a resonant frequency Ω res = Ω 2 0 -ζ 2 2 . Thus, the quantity 2λ can be interpreted as the damping constant ζ of the mechanical system. Therefore, provided that 2λ is much smaller

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

than Ω 0 (i.e. λ 2 /Ω 2 0 1), the model's dynamics is essentially the same as for a linear damped harmonic oscillator with a resonant frequency Ω res = Ω 0 . Yet, while the driving frequency has to be tuned to achieve resonance in the mechanical system, this is not the case for the noisy system that we are considering. In fact, as the system is driven by internal white noise (demographic stochasticity) which covers all frequencies, the resonant frequency of the system will be excited without need of any external tuning.

Other fruitful quantities to measure repetitiveness and (quasi-)periodicity in a fluctuating population are the autocorrelation functions

x A (τ + t)x A (τ ) = x B (t + τ )x B (t) and x A (τ + t)x B (τ ) = -x B (τ + t)x A (τ ) .
The latter are related to the the power spectrum by the Wiener-Khinchin theorem and can be computed by the Fourier transform of P (Ω) [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]), which yields (for t ≥ τ 1)

x A (τ + t)x A (τ ) = 4(1 + 3μ) 3N e -|λ|t |λ| cos (ω 0 t), x A (τ + t)x B (τ ) = 2(1 + 3μ) 3N e -|λ|t |λ| √ 3 sin (ω 0 t) -cos (ω 0 t) . ( 17
)
Following Ref. [START_REF] Nisbet | Modelling fluctuating populations[END_REF]), fluctuations are non-cyclic if the auto-correlation functions decay monotonically to zero, whereas they lead to quasi-cycles if the auto-correlations oscillate. In the latter case, quasi-cycles are "phase-forgetting" if the oscillations of the auto-correlation functions are damped, as in Eqs. ( 17), and "phase-remembering" when their amplitude do not vanish. Following this classification, it is clear from (17) and Figs. 5 and 6 that the quasi-cycles in the region (a) of the parameter space are "phase-forgetting", while those in the region (b) are "phase-remembering".

A analysis similar to that of this section can be carried out in the region (b) of the phase parameter by performing a van Kampen expansion around the limit cycle σ(t). This leads to deal with a stochastic version of non-autonomous linear differential equations (see, e.g., [START_REF] Boland | How limit cycles and quasi-cycles are related in systems with intrinsic noise[END_REF])).

Average escape time

Another intriguing effect of demographic noise is related to the pivotal concepts of attractor and stability.

As discussed in the previous sections there are two different types of stable attractors in the model RPSM:

s * = (1/3, 1/3, 1/3
) is the only attractor of the system when λ ≤ 0, whereas all trajectories in the phase portrait approach the limit cycle σ = (ā(t), b(t)) when λ > 0. This scenario predicting drastically different fate depending on whether λ ≤ 0 or λ > 0 has to be revised when fluctuations are accounted. In this case, within a probabilistic setting, attractors can be regarded as minima of a potential well from which it is always possible to escape after some characteristic time (which can be enormously long, see, e.g., [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]; [START_REF] Risken | The Fokker-Planck Equation[END_REF]; [START_REF] Kubo | Fluctuation and relaxation of macrovariables[END_REF]; [START_REF] Dykman | Large fluctuations and optimal paths in chemical kinetics[END_REF]Volovik et al. (2009))). Stochasticity affects the concept of stability even if the system is initially at the fixed point s * . In fact, due to demographic fluctuations, the species densities deviate from s * and yield quasi-cycles resulting in persistent erratic oscillations of non-vanishing amplitude. It is therefore biologically relevant to assess the robustness of the population composition and understand how noise affects the co-evolution of sub-populations that initially coexist with the same density. Here, to further assess the influence of intrinsic fluctuations when N is large yet finite, we are interested in the average time to escape from s * and reach a specified separating distance from it. In other words, we compute the average time T esc (R) for a trajectory starting at s * to reach a cycle C(R) on which oscillations around s * are of amplitude 2R/ √ 3 (see below). Below, we discuss our results, that are summarized in Fig. 9, and also present an analytical approach -based on the diffusion theory and van Kampen expansion -to compute T esc (R), that is an accurate Typical results are reported in Fig. 9, where one essentially distinguishes two cases. If λ < 0, s * is stable and ln (NT esc (R)) grows linearly with N (for large N ) with a slope depending on R, as illustrated in Fig. 9 (left), which yields to T esc (R) ∝ e C1(R)N /N . On the other hand, when s * is unstable (λ > 0) and N is large, we find that exp(2λT esc (R)) varies linearly with N (slope depending on R), as shown in Fig. 9 (right). This leads to T esc (R) ∝ ln (C 2 (R) N )/(2λ). In the above asymptotic expressions, C 1 (R) and C 2 (R) are monotonic functions of R and comparison with analytical calculations around s * suggests that C 1 ∼ C 2 ∝ R 2 (see below). For the marginal case λ = 0 and large N , numerical results are reported in the inset of Fig. 9 (left) and we find that (ln T esc (R))/R 4 (approximately) exhibits a linear dependence on N (with constant slope). In the case λ = 0, one thus finds T esc (R) ∼ e C3R 4 N , where C 3 is a constant.

In the realm of the diffusion theory and van Kampen expansion, the calculation of the (approximate) average escape time T esc (x) can be regarded as a first-passage problem to an absorbing cycle C(R) starting from a system initially at s * [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]; [START_REF] Redner | A Guide to First-Passage Processes[END_REF]). To obtain an equation for T esc it is useful to consider the socalled backward Kolmogorov equation (BKE), which is the adjoint of Eq. ( 11) and reads (see, e.g., [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]; [START_REF] Risken | The Fokker-Planck Equation[END_REF]; [START_REF] Redner | A Guide to First-Passage Processes[END_REF])):

-∂ t P b (x, t) = L b P b (x, t), with L b (x) = A ij (s * )x j ∂ xi + 1 2 B ij (s * )∂ xi ∂ xj . ( 18
)
As in Section 3.2, to reveal the polar properties of the system in the vicinity of s * , it is natural to perform the (2006)). Thus, in the y-variables, the differential operator of Eq. ( 18) becomes

L b (x) → L b (y) = λ(y A ∂ yA + y B ∂ yB ) + ω 0 (y B ∂ yA -y A ∂ yB ) -1+3μ 6N (∂ 2 yA + ∂ 2 yB )
. To exploit the system's symmetry around s * , we then adopt the polar coordinates (ρ, φ), with y A = ρ cos φ and y B = ρ sin φ. The radial BKE thus reads

∂ t P b (ρ, t) = L b (ρ) P b (ρ, t), with L b (ρ) = λρ + 1 + 3μ 6N 1 ρ ∂ ρ + 1 + 3μ 6N ∂ 2 ρ ( 19
)
To derive the above radial operator, we have assumed that the initial radial symmetry of the probability distribution (initially the system is at s * , without any angular dependence) is approximately preserved by the dynamics. The above approach is certainly valid at linear order around s * and is consistent with the van Kampen linear expansion that has led to (11) and ( 18). With (19), the equation for the average escape time T esc (ρ) at a distance ρ from s * is given by [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]; [START_REF] Risken | The Fokker-Planck Equation[END_REF])

-1 = L b (ρ) T esc (ρ). ( 20
)
In the framework of diffusion and van Kampen approximations, we now analytically compute the mean escape time T esc (R) for ρ to attain the value R starting from s * (i.e. initially ρ = 0), which corresponds to the mean time necessary to reach the cycle C(R). Our analytical treatment is valid in the linear regime around s * , where φ ω 0 , and the amplitude of the oscillations on C(R) is 2R/ √ 3. In our theoretical setting, the average escape time T esc (R) is obtained by solving Eq. ( 20) subject to a reflecting and an absorbing boundary conditions at ρ = 0 and ρ = R, respectively [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]). When λ = 0, the solution to this problem can be expressed in terms of the function Ψ(z) = z exp 3Nλ 1+3μ z 2 :

T esc (R) = 6N 1 + 3μ R 0 dy Ψ(y) y 0 dzΨ(z) = 1 2λ -3NR 2 λ 1+3μ 0 du u (1 -e u ) . ( 21
)
The solid curves in Fig. 9 show the comparison between the prediction (21) and numerical results. We notice that for small values of R, e.g., R = 0.10 -0.15, there is an excellent agreement between ( 21) and the results of numerical simulations (both in the cases λ < 0 and λ > 0). For larger values of R (e.g. for R ≥ 0.17), nonlinearities and non-vanishing angular dependence cause some systematic deviations from the results of numerical computations (discrepancies of ≈ 5% -10% in the results of Fig. 9 for R = 0.17), but (21) still provides the correct (qualitative) behavior and a reasonable approximation of T esc (R). In the asymptotic limit of large population size N , two distinct behaviors can be obtained from Eq. ( 21):

• When λ < 0, s * is stable and the main contribution to ( 21) is given by T esc (R)

1 2|λ| Ei 3|λ|R 2 N 1+3μ
, where

Ei(x) ≡

x -∞ dt (e t /t) is the exponential integral (Abramowitz and Stegun (1965)). From the properties of this function, we infer the following asymptotic behavior (for |λ|NR 2 1):

T esc (R) 1 + 3μ 6(λR) 2 N exp 3|λ|R 2 N 1 + 3μ (22)
This result predicts that the mean escape time increases dramatically with the total number of individuals N and, when |λ|NR 2 1, is proportional to the exponential of the population size. The dashed lines in Fig. 9 (left) illustrate that ( 22) is a very good approximation of T esc (R) in the linear regime around s * , e.g.

for R = 0.10 -0.15, when N 13 . In particular, we notice that (22) coincides with the aforementioned asymptotic expression of T esc (R) with C 1 = 3|λ|R 2 /(1 + 3μ). For R = 0.17, when N is large, nonlinear effects are important and in Fig. 9 the exponent of ( 22) overestimates that of T esc (R) by ≈ 5%.

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

• When λ > 0, demographic fluctuations always cause the departure from s * (unstable), and all trajectories in the phase portrait are attracted by the limit cycle σ. In this case, from the properties of the exponential integral (Abramowitz and Stegun (1965)), for λN R 2 1, the main contribution to ( 21) is

T esc (R) 1 2λ ln 3λR 2 N 1 + 3μ + γ EM , ( 23
)
where γ EM = 0.57721... is Euler-Mascheroni's constant. This result predicts that, when λN R 2 1, the mean escape time grows logarithmically with the population size N , i.e. T esc (R) → (ln N )/2λ, and the value of R contributes to subdominant corrections. Here also, the dashed lines in Fig. 9 (right), almost indistinguishable from solid curves, agree very well with the numerically computed T esc (R) around s * (i.e.

for R = 0.10 -0.15). We notice that, when C 2 = 3λR 2 /(1 + 3μ), the prediction (23) coincides with the asymptotic expression of T esc (R) inferred from Fig. 9 (right).

Besides the above cases, a special situation arises when λ = 0, μ > 0 and s * is a center. In fact, in Sec. 3.2, we have seen that in this case s * is an attractor, but the dynamics in its vicinity is very slow. To account for such a slow dynamics it is not sufficient to consider (11, 18) -obtained by linearization about s * -but nonlinearities have to be accounted. Generalizing the above approach and keeping the cubic terms in the drift contribution [START_REF] Cremer | Anomalous finite-size effects in the Battle of the Sexes[END_REF]), leads to (20) with the following differential operator [instead of

(19)]: L b (ρ) = {λ + βρ 2 }ρ + 1+3μ 6N 1 ρ ∂ ρ + 1+3μ 6N ∂ 2 ρ .
When λ = 0, proceeding as above, one still finds an exponential asymptotic behavior (μR 4 N 1): T esc (R) ∼ N -1/2 exp(27μR 4 N/(1 + 3μ)). This estimate, that is consistent with the numerical results reported in Fig. 9 (left, inset), implies that (at given R and N ) the mean escape time T esc around s * is significantly shorter in the marginal case than in the case λ < 0.

Summary and Conclusion

This work has been motivated by the importance of further understanding the mechanisms at the origin of persistent and erratic oscillatory dynamics in population biology, which is frequently observed but whose theoretical origin is often debated. Furthermore, this work concerns the effects of mutations in rock-paper-scissors (RPS) games, which are regarded as paradigmatic models to describe the co-evolutionary dynamics of systems with codominant interactions. It has recently been proposed that codominance can help promote the maintenance of biodiversity, and it is therefore of biological relevance to study generalizations of the RPS model to understand which ingredients favor the long-term coexistence of the sub-populations.

Here, we have investigated the oscillatory dynamics of the generic rock-paper-scissors games with mutations in a well-mixed (homogeneous) population of N individuals. Our study has been carried out in the mean-field limit (N → ∞) and in the presence of demographic noise (N large but finite). In addition to the regions of the parameter space associated with a stable interior fixed point s * and heteroclinic cycles, the possibility for the individuals to switch from one strategy (species) to another with a small transition rate yields a new oscillatory behavior associated with a limit cycle and resulting from a Hopf bifurcation. In the mean-field limit (N → ∞), we have recast the system's rate equations in a normal form and studied the main properties of the ensuing limit cycle. When the population size is finite, demographic noise has to be taken into account and the dynamics becomes genuinely stochastic and nonlinear. The latter situation has been described in terms of an individual-based formulation, with the stochastic dynamics implemented according to a Moran process. The influence of intrinsic noise on the dynamics has been analytically assessed within a diffusion theory approximation complemented by numerical simulations (using the Gillespie algorithm). We have shown that demographic noise transforms damped oscillations into ("phase-forgetting") quasi-cycles and perturbs the amplitude and the phase of limit cycles (that become "phase-remembering" quasi-cycles). As also observed in other systems (e.g., [START_REF] Bartlett | Measles periodicity and community size[END_REF][START_REF] Bartlett | The critical community size for measles in the United States[END_REF]; [START_REF] Nisbet | Modelling fluctuating populations[END_REF]; [START_REF] Mckane | Predator-prey cycles from resonant amplification of demographic stochasticity[END_REF]; [START_REF] Boland | How limit cycles and quasi-cycles are related in systems with intrinsic noise[END_REF])), the effect of stochasticity is particularly striking in the region of the parameter space where s * is stable. There, we have found sustained oscillations driven by demographic stochasticity. In fact, while fluctuations are of order N -1/2 , for small mutation rate, the amplitude of the oscillations is amplified by a resonance amplification caused by internal (white) noise. We have shown that this phenomenon translates into an isolated and wellmarked peak in the power-spectrum and have also computed the autocorrelation functions of the system. To further assess the robustness of the co-evolutionary dynamics against noise, we have computed numerically and analytically -using the diffusion approximation and the mapping onto a first-passage problem -the average escape time necessary to reach a cycle on which the oscillations attain a given amplitude. We have therefore shown that such a mean exit time grows either exponentially or logarithmically with the system size, depending on whether the interior fixed point is deterministically stable or unstable, respectively.

In summary, in the presence of mutations the RPS model yields sustained oscillations in every region of the parameter space. The latter are either caused by demographic noise and/or result from a Hopf bifurcation.

Therefore, mutations in the RPS model ensure that all species co-evolve and oscillate in time without ever going extinct. One can therefore regard mutations as a mechanism favoring the oscillatory dynamics in communities characterized by cyclic co-dominance.

In future research, it would be interesting to consider a spatially-extended version of this model and study how the population will self-organize and which kinds of spatial-temporal patterns would emerge. In particular, as it has recently been shown that mobility can affect co-evolution in RPS models by mediating between meanfield and stochastic dynamics (Reichenbach et al. (2007a(Reichenbach et al. ( ,b, 2008b))), it would be relevant to investigate variants of the model RPSM where individuals are allowed to move and interact in space.

the

  Moran process (Moran (1958);[START_REF] Nowak | Evolutionary Dynamics[END_REF]). Recently, individual-based modeling has been extensively used to study the stochastic RPS dynamics. In fact, a large body of recent research has been concerned with spatially-extended systems and the influence of the species' movement on the emerging noisy patterns (Tainaka

  a and π = (1 -)(ab + bc + ac). The latter are studied in the next sections.

Figure 1 :

 1 Figure 1: (Color online). Typical plots of the time-dependent density (relative abundance) a(t) of species A, obtained from Eqs. (2) for = 1.225 and different values of μ (with a(0) = 0.25, b(0) = 0.40).Here, μ c = 0.0125 (see text). Left: In the absence of mutation rate (μ = 0, λ = 0.0375), the fixed point s * is unstable and the density (brown, dashed) jumps from 0 to 1. For very low values of the mutation rate, here μ = 0.001 (λ = 0.0345), s * is also unstable and a(t) oscillates regularly (magenta, solid) with an amplitude that approaches the extreme values 0 and 1, but with a much shorter period than in the case μ = 0. Right: For μ = 0.01 (λ = 0.0075), the fixed point s * is still unstable and a(t) oscillates regularly about a * = 1/3 (red, solid), with a period T ≈ 10. When μ = 0.02 > μ c (λ = -0.0225), the densities exhibit exponentially damped oscillations and converge towards the fixed point value 1/3 (blue, dashed).

Figure 2 :

 2 Figure 2: (Color online). Bifurcation diagram in the parameter space of the (deterministic) model RPSM. The regions (a) and (b) are separated by the critical mutation rate μ c = ( -1)/18 > 0. In region (a), μ > μ c for > 1 and μ ≥ 0 for ≤ 1, s * is a stable focus. Below the line μ c ( ), in region (b)where 0 < μ < μ c , the system undergoes a Hopf bifurcation (at μ = μ c ) and the dynamics is characterized by a limit cycle and stable oscillations. In region (c), where μ = 0 (no mutations) and > 1 (dashed line), the REs (2) yield heteroclinic cycles along with periodic oscillations unstable against demographic fluctuations.

Figure 3 :

 3 Figure 3: (Color online). Plot of r ∞ = λ/|β| > 0 for μ > μ c as function of and μ according to (4) and (7). We notice that r ∞ increases when is raised and μ is lowered.

Figure 4 :

 4 Figure 4: (Color online). With xA (t) = ā(t) -1 3 and xB (t) = b(t) -1 3 , comparison of the limit cycle obtained from the numerical solution of (2) (red/thick, here 550 ≤ t ≤ 1000) with the predictions of (6) [black/thin]. The parameters and initial conditions are = 1.2, μ = 0.01, x A (0) = 0.025, x B (0) = 0.040, here λ = 1/300.

Figure 5 :

 5 Figure5: (Color online). Erratic oscillations and stochastic dynamics of the model RPSM in the region (a) of the parameter space, with = 0.5, μ = 0.01 and N = 300 (initially N A = 135 and N B = 66). In single realizations, the time-dependent densities (relative abundances) of species A (thin, red/dark grey) and B (thin, cyan/light grey) are characterized by persistent stochastic oscillations resulting from a resonance amplification occurring at frequency Ω * ω 0 0.43 (see text). When they are sample-averaged (here over 200 replicates), the mean densities (solid and dashed thick curves for species A and B, respectively) display damped oscillations approaching the value 1/3, in agreement with (2).

  Figs. 5 and 6 show that demographic noise affects considerably each single realization in the region (a) of the parameter space, where the predictions of the REs are replaced by oscillations of non-vanishing amplitude, while it perturbs the cyclic behavior in region (b). Below, we show that the "quasi-cycles" arising in region (a) of Fig.2at low mutation rate result from a resonance amplification mechanism[START_REF] Mckane | Predator-prey cycles from resonant amplification of demographic stochasticity[END_REF]

Figure 6 :

 6 Figure 6: (Color online). Stochastic dynamics of the model RPSM in the region (b) of the parameter space with = 1.25, μ = 0.008 and N = 600 (initially N A = N B = 240). We notice that the averaged (over 250 replicates) time-dependent densities of species A (thick, solid) and B (thick, dashed) are characterized by oscillations of non-vanishing amplitude. The erratic trajectories in red/dark grey and cyan/light grey are single realizations (no sample averaging) of the density of species A and B, respectively.

Figure 7 :

 7 Figure 7: (Color online). Phase portrait of the model RPSM in the presence of stochastic fluctuations for = 1.25, N = 600, with μ = 0.008 (thick, red/dark grey) and μ = 0.05 (thin, cyan/light grey), initially N A = N B = 240. The trajectories have been obtained from a sample-average over 250 replicates.When μ = 0.008, the interior fixed point s * is unstable and each trajectory in the phase portrait forms a "phaseremembering quasi-cycle"of erratic non-vanishing amplitude and phase orbiting around σ. When μ = 0.05, the flow spirals towards the stable interior fixed point s * and erratically wanders in its vicinity forming a "phase-forgetting quasi-cycle" (see

Figure 8 :

 8 Figure 8: (Color online). Power spectrum P (Ω) = | x(Ω)| 2 of the model RPSM with = 0.5, μ = 0.01 and N = 300 obtained from numerical simulations and by sample averaging over 200 realizations. We notice an isolated peak at characteristic frequency Ω * 0.43, in full agreement with the analytical prediction (15).

  resonant amplification of demographic (intrinsic) fluctuations (McKane and Newman (2005)).

Figure 9 :

 9 Figure 9: (Color online). Average escape time T esc (R) to reach the cycles C(R) starting from s * , for R = 0.10 (♦), 0.12 ( ), 0.15 (•), 0.17 (×) and large population size N . Left: Plot of ln (NT esc (R)) versus N in the case λ < 0, with = 0.5 and μ = 0.01, and sample-average over 100 replicates. We find a linear dependence that yields T esc (R) ∝ e C1(R) N /N . Numerical results (filled symbols) are also compared to the predictions of Eqs. (21) (solid) and (22) (dashed), yielding C 1 (R) ∝ R 2 . Right: Plot of exp(2λT esc (R)) versus N in the case λ > 0, with = 1.5 and μ = 0.01 and sample-average over 400 replicates. The observed linear dependence yields T esc (R) ∝ ln (C 2 (R) N )/(2λ). Numerical results (filled symbols) are also compared to the predictions of Eqs. (21) (solid) and (23) (dashed lines, almost indistinguishable from the solid lines), yielding C 2 (R) ∝ R 2 . Inset: Marginal case λ = 0 with parameters = 1.18 and μ = 0.01, and sample-average over 100 replicates. Plot of R -4 ln (T esc (R)) versus N for values of R = 0.30 ( ), 0.35 ( ), 0.40 (∇), 0.45 (+). For large population size, the (approximate) linear dependence R -4 ln (T esc (R)) ∝ 0.23N (dashed line as a guide for the eyes) yields T esc (R) ∼ e 0.23R 4 N (see text).

  approximation of T esc (R) around s * . Using the Gillespie algorithm, we have computed numerically the average escape time T esc (R) necessary to reach a cycle C(R) starting from the interior fixed point a * = b * = c * = 1/3 (same initial density of each species). To exploit the symmetry of the model (see Sec. 3.2), we work with the variables y = (y A , y B ) = Sx and consider cycles C(R) of parametric equation y A = R cos φ, y B = R sin φ. In the original x-variables, the equation of C(R) is x A = R cos φsin φ/ √ 3 , x B = (2R/ √ 3) sin φ and the amplitude of the oscillations is 2R/ √ 3. In our numerical computations, we have considered systems of population size N , with N = 300 -4000, and values of R ranging from 0.10 to 0.45. Each simulation has been sample-averaged over 100 to 400 replicates.As we are mainly interested in the linear regime around s * , where analytical calculations are amenable (see particular attention has been dedicated to values of R ≈ 0.10 -0.20 (i.e. amplitude ≈ 0.11 -0.23).

  linear transformation x → y = Sx and work in the y-variables. Under this transformation, one has P b (x, t) → P b (y, t), A(s * ) → SA(s * )S -1 =

In fact, the number of individuals of species A is related to a by N A = aN + O( √ N ) and similarly with the other species.

E.g., for = 0.5 and μ = 0.01, one finds: ω0 0.4330, whereas Ω * 0.4328 and Ω0 0.4476.

In Fig.9(left), 3|λ|R 2 N/(1 + 3μ) lies between ≈ 1.0 (for R = 0.1 and N = 300) and ≈ 9.5 (for R = 0.17 and N = 1000).

Acknowledgments

This work was initiated with the partial support of the Swiss National Science Foundation through the Advanced Fellowship PA002-119487.

A c c e p t e d m

a n u s c r i p t Szabó, G., Szolnoki, A., 2002


Joseph Hardwicke

Bing Song

Ryan Moseley

David W Thomas

INVESTIGATION OF THE POTENTIAL OF POLYMER THERAPEUTICS IN CORNEAL RE-EPITHELIALISATION

Keywords: cornea, wound healing, animal study, pharmacology

We describe the first use of a bioresponsive polymer therapeutic agent in the promotion of corneal re-epithelialisation after injury in an ex vivo whole-eye organ culture model. A polymer-protein conjugate consisting of dextrin and recombinant human epidermal growth factor (rhEGF) was synthesized and applied as a single dose to a 2mm ex vivo corneal ulcer, in culture. Enhanced wound healing was observed in response to dextrin-rhEGF, when exposed to α-amylase, compared to controls. This highlights the potential for polymer therapeutics to provide a platform for bioresponsive drug/protein delivery in the field of ophthalmology.

INNOVATION:

The descriptor "polymer therapeutic" is an umbrella term used to describe polymeric drugs, polymer-drug conjugates, polyplexes, polymer-protein conjugates and polymeric micelles, to which a drug can be bound [START_REF] Duncan | The dawning era of polymer therapeutics[END_REF] . In the context of polymer-protein conjugates, the stability and solubility of proteins can be enhanced, and their biological half-life and immunogenicity reduced, via polymer conjugation [START_REF] Veronese | Introduction and overview of peptide and protein pegylation[END_REF][START_REF] Veronese | PEGylation, successful approach to drug delivery[END_REF] .

To this date, although polymer therapeutics have been licensed for use in humans (Table 1), they have not been investigated widely in the field of corneal wound healing.

We have been developing a novel controlled release formulation of recombinant human EGF (rhEGF), conjugated to the degradable polysaccharide, dextrin [START_REF] Hardwicke | Dextrin-rhEGF conjugates as bioresponsive nanomedicines for wound repair[END_REF] . This provides a controlled delivery platform for the release of rhEGF in response to endogenous α-amylase, which is present in the tear film [START_REF] Van Haeringen | Amylase in human tear fluid: Origin and characteristics, compared with salivary and urinary amylases[END_REF] . When conjugated to the dextrin polymer, the bioactivity of the growth factor is decreased, but an increase in protein stability is afforded. Bioactivity can be fully restored after polymer degradation by α-amylase, and the dextrin breakdown products, maltose and isomaltose are further metabolised and/or excreted [START_REF] Peers | Icodextrin provides long dwell peritoneal dialysis and maintenance of intraperitoneal volume[END_REF] . This Polymer-masking UnMasking Protein Therapy hypothesis PUMPT hypothesis [START_REF] Duncan | Polymer Masked-Unmasked Protein Therapy (PUMPT) 1. Bioresponsive dextrin-trypsin and -MSH conjugates designed for α-amylase activation[END_REF] , has been shown to stimulate dermal keratinocyte proliferation in vitro, and enhances both the half-life of rhEGF and the stability in response to wound proteolysis [START_REF] Hardwicke | Dextrin-rhEGF conjugates as bioresponsive nanomedicines for wound repair[END_REF] .

In animal studies, EGF stimulated proliferation and DNA synthesis of corneal epithelium and increased the epithelial healing rate [START_REF] Frati | Selective binding of the epidermal growth factor and its specific effects on the epithelial cells of the cornea[END_REF][START_REF] Savage | Proliferation of corneal epithelium induced by epidermal growth factor[END_REF][START_REF] Ho | Kinetics of corneal epithelial regeneration and epidermal growth factor[END_REF][START_REF] Watanabe | Stimulatory effects of fibronectin and EGF on migration of corneal epithelial cells[END_REF][START_REF] Singh | Epidermal growth factor in alkali burned corneal epithelial wound healing[END_REF][START_REF] Gönül | Effect of EGF dosage forms on alkali burned corneal wound healing of mice[END_REF] . The initial hope of EGF to progress to a licensed product for clinical use in the early 1990's [START_REF] Caporossi | Epidermal growth factor in topical treatment following epikeratoplasty[END_REF][START_REF] Scardovi | Epidermal growth factor in the topical treatment of traumatic corneal ulcers[END_REF][START_REF] Cellini | Epidermal growth factor in the topical treatment of herpetic corneal ulcers[END_REF] , did not come to fruition. The lack of significant effects from EGF, in humans, applied post penetrating keratoplasty [START_REF] Dellaert | Influence of topical human epidermal growth factor on postkeratoplasty re-epithelialisation[END_REF] , did not match the previous findings from animal studies. The aforementioned clinical studies have involved the repeated administration of EGF, due to the short half-life and the action of the tear film [START_REF] Chan | Pharmacokinetic study of recombinant human epidermal growth factor in the anterior eye[END_REF] .

METHODS:

We have evaluated the efficacy of single dose dextrin-rhEGF in an ex vivo organ culture model of corneal wound healing [START_REF] Ramaesh | Increased apoptosis and abnormal wound healing responses in the heterozygous Pax6 +/mouse cornea[END_REF][START_REF] Reid | Wound healing in rat cornea: The role of electric currents[END_REF] . The synthesis and characterization of the conjugate was reproducible (Figure 1) [START_REF] Hardwicke | Dextrin-rhEGF conjugates as bioresponsive nanomedicines for wound repair[END_REF] . Two-month old male Wistar Han rats (250-300g; Charles River Laboratories, Margate, UK), were sacrificed by CO 2 asphyxiation, confirmed by cervical dislocation. Wounding was performed within 1 hour of sacrifice. The globe was irrigated with Balanced Salt Solution (BSS Plus ® ; Alcon Laboratories, Hemel Hemstead, UK) containing 0.1% EDTA. A full-thickness 2 mm corneal abrasion was produced by excision. The globe was then enucleated. Tissue culture medium (DMEM supplemented with penicillin G, 100U/mL; streptomycin sulphate, 100µg/mL; amphotericin B, 0.25µg/mL; Lglutamine, 2mM (Invitrogen; Paisley, UK): serum-free media; SFM) was supplemented with study compounds, as shown in Table 2, and added to the globes (n = 4). Specimens were maintained at 37 °C / 5% CO 2 , for 64 h. Wound healing was assessed by digital wound photography and fluorescein staining. After 64h, globes were fixed in 10% formal saline, embedded in wax, sectioned and stained with H&E.

Statistical analyses were undertaken using GraphPad Prism ® , version 4.00 (GraphPad Software, San Diego, USA). Data were compared using a one-way ANOVA, with a Bonferroni post-test, for group analysis. Results were expressed as a mean and standard deviation (S.D.). Statistical significance was considered at a probability of P < 0.05.

RESULTS:

This ex vivo model of corneal wound healing was sensitive to stimulation by foetal calf serum (FCS) and by dextrin-rhEGF. In the conjugated form, the rhEGF promoted a significant increase in the rate of corneal wound healing, with wound closure by 48 hours, in the presence of human physiological levels of α-amylase [23][24] compared to free rhEGF. With the addition of exogenous α-amylase, at physiological concentrations, rhEGF was released, or "unmasked" from the polymer, thereby restoring bioactivity (Figures 2 and3).

There was a significant improvement in wound healing (P < 0.006), in response to serum (SCM; Group 2), compared to the serum-free (SFM; Group 1) and α-amylase (Group 3) controls, at 24 h (Figure 4a). All wounds in the SCM group (Group 2) had fully re-epithelialised by 48 h, whilst the average wound areas in the serum-free controls (Groups 1 and 3), were 24.5 ± 8.1% and 31.4 ± 5%, respectively, of that at time 0 h (P < 0.001). In the unconjugated rhEGF Groups (Groups 4 and 5), there were no significant effects on wound re-epithelialisation (P > 0.05), at either concentration, and neither resulted in total wound re-epithelialisation. Wound areas, at 64 h, were 5.6 ± 11.3% and 23.4 ± 21.4%, of the total wound area at time 0 h, for the 10 μg/mL rhEGF, and 1 μg/mL rhEGF Groups, respectively (mean ± S.D.) (Figure 4b).

In the dextrin-rhEGF conjugate Study Groups (Groups 6 to 8), there was failed re-epithelialisation, in the absence of physiological concentrations of α-amylase (Group 6; "masked"), with an abrasion of 23.4 ± 5.5% of that at time 0, persisting at 64 h. This was not significantly different to the SFM control (Group 1) (P > 0.05).

Upon the addition of α-amylase (93 i.u./L; "unmasking"), there was a significant decrease in the wound area, at 24 h, in both the 1 μg/mL rhEGF equivalent (Group 7) and 10 μg/mL rhEGF equivalent (Group 8) Groups (P < 0.006). The dextrin-rhEGF conjugate was significantly more effective, at decreasing the wound area, at 24 h, at the lower dose (1 μg/mL rhEGF equivalent), than at the higher dose (10 μg/mL rhEGF equivalent) (P = 0.027) (Figure 4c).

DISCUSSION:

The increased efficacy of the dextrin-rhEGF conjugate, in enhancing corneal wound healing, may be due to the sustained release phenomena [START_REF] Hardwicke | Dextrin-rhEGF conjugates as bioresponsive nanomedicines for wound repair[END_REF] . In the ex vivo model described herein, the rhEGF was applied as a single dose at the outset of the experiment, and the free rhEGF may have been degraded, due to the short half-life of EGF [START_REF] Chan | Pharmacokinetic study of recombinant human epidermal growth factor in the anterior eye[END_REF] . Although initial wound closure rates are similar in the first 24 h between rhEGF and dextrin-rhEGF (+α-amylase), this apparent difference tailed off as time progressed. In the ex vivo organ culture model, the pharmacokinetics of EGF will differ, but the sustained release of rhEGF, from the α-amylase-mediated, degradation of the dextrin-rhEGF conjugate, may account for the enhanced wound healing in the dextrin-rhEGF study groups. The ex vivo model of corneal wound healing, has been shown to be reproducible, and that it is a reliable model for further analysis of polymer-protein conjugates. This is the first application of polymer therapeutics specifically designed for tissue regeneration, in such a model, and these preliminary investigations highlight the potential for the application of these therapeutics to enhance corneal wound healing. This application of a specific endogenous trigger to promote the controlled release of a bioactive peptide, may be applied to the healing and/or regeneration of other tissues in the future. Polymer therapeutics that are licensed in the UK. PEG = polyethylene glycol, GCSF = granulocyte colony stimulating factor, IFN = interferon.

Table 2

Culture medium supplementation. PBS = phosphate buffered saline, rhEGF = recombinant human epidermal growth factor.

The PUMPT hypothesis. Succinoylated dextrin is conjugated to rhEGF using EDC and sulfo-NHS cross-linkers, to produce a dextrin-rhEGF polymer-protein conjugate.

On exposure to α-amylase (present in human tear fluid), bioactive rhEGF is released from the polymer in a controlled release fashion. rhEGF would induce corneal epithelial cell proliferation and migration, whilst the polymer would be metabolised to the maltose and isomaltose disaccharide units. RT = room temperature.

Figure 4

Panel (a) shows the corneal wound healing (% reduction in wound area compared to 0 h), over 64 h, in response to FCS (serum-containing media) and serum-free controls.

There was a significant decrease (*P <0.006) in wound area, in response to FCS, when compared to serum-free controls. Panel (b) shows the corneal wound area reduction, in response to free rhEGF. The initial rate of wound area reduction (over the first 24 h), was not significantly different (P > 0.05) from FCS, but the effect was not sustained over 64 h. Panel (c) shows the corneal wound area reduction, in response to dextrin-rhEGF. In the "masked" (α-amylase-free) state, the bioactivity of dextrin-rhEGF was reduced to that of serum-free controls (P > 0.05), but upon the addition of α-amylase, at physiological concentrations, bioactivity was significantly restored (**P < 0.006). Dextrin-rhEGF was most effective at the 1 μg/mL concentration (***P = 0.027), compared to the conjugate at 10 μg/mL concentrations.

Figure 2 Typical

 2 Figure 2

Figure 3 Panel

 3 Figure 3

Table 1

 1

	Name	Manufacturer	Polymer	Protein	Indication
	Neulasta™	Amgen	PEG	GCSF	Chemotherapy-
	(Pegfilgrastim)				induced
					neutropaenia
	Pegasys™	Roche	PEG	α-IFN-2a	Hepatitis C
	(Peginterferon alfa)
	ViraferonPeg™	Schering-Plough	PEG	α-IFN-2b	Hepatitis C
	(Peginterferon alfa)
	Somavert™	Pfizer	PEG	Human	Acromegaly
	(Pegvisomant)			growth
				hormone
				antagonist

TABLE LEGENDS Table 1

 LEGENDS1

ACKNOWLEDGEMENTS: J.H. was funded by a Fellowship from The Healing Foundation (Registered Charity number 1078666) and the Welsh Office for Research and Development (WORD) and the authors would like to acknowledge support from EPSRC platform grant No. EP/C 013220/1.

COMPETING INTEREST: None declared.

COPYRIGHT STATEMENT:

The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive license (or non-exclusive for government employees) on a worldwide basis to the BMJ Publishing Group Ltd and its Licensees to permit this article (if accepted) to be published in British Journal of Ophthalmology and any other BMJPGL products to exploit all subsidiary rights, as set out in the license.