Edward Malone

Chris Elliott

David Glenn Kennedy

David Savage

Liam Regan

Glenn Kennedy

Surveillance study of a

Keywords: Veterinary drug residues -anabolic steroids Bovine tissue, Mass Spectrometry, Steroids, Imported, Survey

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In recent years there has been an increase in global concern relating to the adverse effects that may result from exposure to chemicals that could interfere with the endocrine system. A number of classes of compound have been highlighted in this regard, they include dioxins, phyto-oestrogens, Bisphenol A, hormones, etc. of these hormones continue to be in general the most biologically active. Hormonally active growth promoters have been shown to be very effective for promotion of muscle growth in farm animals; particularly bovine ruminants, veal calves and pigs [START_REF] Lu | Enviromental, quality and safety[END_REF], WHO 1982, Meissonnier et al. 1983[START_REF] Berende | Domestication, conservation and use of animal resources[END_REF]. Within the EU, treatment of animals intended for human consumption with substances exhibiting thyreostatic, androgenic, oestrogenic or gestagenic effects is a regulatory offence (Council Directive EC 96/22 1996). The use of these substances is controlled through national residue control plans in each EU member state (Council Directive EC 96/23 1996). The EU inspection programs focus on matrices that are more suitable for the detection of prohibited substances, especially if the animals are still on farm, these matrices urine, faeces or hair. With regard to slaughtered animals, the more favoured matrices are bile, retina, fat, blood, liver and kidney. Rarely is muscle tissue sampled for the detection of steroidal type growth promoters. Muscle is usually only tested in the case of import controls or in monitoring programs of meat sampled in butcher shops or supermarkets. As a result, data on steroid concentrations in muscle meat samples is rare. Data on levels of natural hormones in meat is even rarer because of the absence of "legal natural levels" for these hormones in compliance testing. A restricted controlled use of some specified anabolic steroids is permitted in the USA, Canada, Australia, New Zealand and in some countries in South America, Asia and Africa. The EU imports large amounts of beef each year from these third countries, 524,000 tonnes in 2005 alone (Eurostat 2005). However when beef is imported in to the EU, it should meet the same standards with regards the absence of prohibited substances as prevail within the EU i.e. absence of illegal growth promoters. As the number of growth promoters likely to be abused in cattle fattening increases the use of multi-residue methods that cover a number of classes of compounds becomes increasingly important. Methods such as this are employed in the testing of domestic samples of urine as part of the national control plan each EU state must carry out. These methods usually monitor for a large amount of analytes for example 22 compounds [START_REF] Hewitt | Screening and confirmatory strategies for the surveillance of anabolic steroid abuse within Northern Ireland[END_REF]. Very sensitive analytical techniques for the detection of these substances and their metabolites in tissue need be adopted as the levels that accumulate in tissue are lower than those found in other matrices (Marchard et al. 2000). Liquid chromatography tandem mass spectrometry (LC-MS/MS) has been the determination method of choice for the majority of articles published recently concerned with the determination of steroidal growth promoters in bovine tissue (Xu et al. 2006, Chursch et al 2008[START_REF] Costain | Identification of hormone esters in injection site in muscle tissues by LC/MS/MS[END_REF][START_REF] Malone | Development of a rapid method for the analysis of synthetic growth promoters in bovine muscle using liquid chromatography tandem mass spectrometry[END_REF], Schmidt et al 2009[START_REF] Kaklamanos | Determination of anabolic steroids in muscle tissue by liquid chromatography-tandem mass spectrometry[END_REF], Yang et al 2009). Various types of extraction and purification protocols have been adopted in these procedures, the majority of which use solid phase extraction (SPE) for sample purification; in some cases a dual SPE step is adopted. A wide range of SPE sorbents have been utilised for this task; Prior to extraction and purification a number of methods employ a deconjugation step in order to release any glucuronide or sulphate forms of the compounds to their free form. This step is required for urine samples as a considerable percentage of these compounds are excreted in the glucuronide and sulphate forms in urine (Le Bizec et al 1992). However in other matrices such as muscle, evidence suggests that this step may not be necessary, investigations done by a number of groups found that when analysing muscle samples no significant difference in concentrations of analyte observed could be found between samples which underwent enzymatic hydrolysis and those which did not [START_REF] Hartmann | Simultaneous determination of anabolic and catabolic steroid hormones in meat by gas chromatography-mass spectrometry[END_REF][START_REF] Marchand | Ultra trace detection of a wide range of anabolic steroids in meat by gas chromatography coupled to mass spectrometry[END_REF], Yang et al. 2009).

The main aims of the work in this study were to develop and validate a rapid confirmatory method for a wide number of steroidal type growth promoters in bovine muscle and apply this method to the analysis of a number of samples of bovine muscle imported into Ireland from third countries. Both natural and synthetic compounds were investigated. However the presence of any naturally occurring compounds such as progesterone, testosterone etc does not necessarily imply misuse.

Experimental

Materials and reagents

Water, methanol, acetonitrile (HiPerSolv grade) were obtained from BDH (Merck, Poole, Dorset, UK). Ammonium acetate and sodium sulphate (ACS grade) were obtained from Sigma (St. Louis, MO, USA). Dexamethasone from QMX Laboratories (Thaxted, UK).

Primary stock standard solutions of each analyte were prepared in methanol at a concentration of 100 µg mL -1 . Intermediate single standards solutions each analyte were prepared in methanol at a concentration of 10 µg mL -1 . A mixed standard fortification solution which included each analyte was prepared in methanol at a concentration of 100 ng mL -1 . A mixed standard fortification solution including each internal standard was prepared in methanol at a concentration of 500 ng mL -1 . All standards were stored at -20 o C. The injection solvent consisted of 0.5mM Ammonium Acetate in methanol:water (30:70, v/v).

LC conditions

The LC system was a Shimadzu UFLC-XR equipped with a LC-20AD-XR Binary pump, SIL-20AD-XR autosampler and a CTO-20A column oven (Shimadzu, Dublin, Ireland). The compounds were chromatographed on a 1.8 µm Agilent Eclipse Plus C18 column (2.1 × 50 mm) (Agilent, Dublin, Ireland) and the column temperature was maintained at 55 o C. A gradient LC system using 0.5 mM ammonium acetate: methanol (70:30, v/v Mobile Phase A) and ammonium acetate: methanol (5:95, v/v Mobile phase B) at a flow of 0.5 mL min -1 was operated. The gradient profile began at 90% A held for 3 minutes before changing to 50% A 2 minutes later, held for 1 minute and changing to 15% A, 1.5 minutes later and held for 1.5 minutes before returning to 90% A after 1.5 minutes and left to equilibrate for 2.5 minutes. The total runtime of the method was 13 minutes.

MS/MS parameters

The mass spectrometer system used was a 5500 triple qaudrupole instrument (Applied Biosystems/ MDS Sciex, Canada).The analysis was performed using both positive and negative ion electrospray (ESI) and the MS was operated in multiple reaction monitoring (MRM) mode. Two transitions per analyte were monitored whereas only one transition was monitored for each deuterated internal standard. The precursor product transitions for each analyte as well as their corresponding collision energies are shown in Table 1. The MS/MS detector conditions were as follows: Ion spray voltage 1500V; Ion Source Gas 1 50 psi; Ion Source Gas 2 60 psi; Temperature 650 o C; Curtain Gas 15 psi; Collision Gas Pressure 8 (arbitrary units); Entrance Potential 10V. Data acquisition and integration were performed using Analyst software version 1.5 (Applied Biosystems/ MDS Sciex, Canada).

Bovine tissue samples

Bovine tissue samples were obtained and stored at -20 °C. Samples of this tissue were analysed and those found to contain no detectable residues of the analytes of interest were used as negative controls to carry out development and validation. Tissue samples (approximately 50 g) were sliced into very fine pieces and mixed together. 3 grams of the resulting mixture is weighed into Foodcon7 TM polypropylene tubes. 30 µL of the mixed internal standard was added at this time, corresponding to a concentration of 5 µg kg -1 . Samples used to calculate the precision and recovery of the method were fortified at levels corresponding to 1, 1.5 and 2 µg kg -1 for each analyte by adding 30, 45 and 60 µL portions of the 100 ng mL -1 mixed standard solution. After fortification, samples were held for 15 min prior to extraction. Acetonitrile (10 mL) was added to the samples followed by sodium sulphate 2.5 grams. The samples were then vigorously shaken (30 seconds) followed by homogenisation using a MediFASTH 2 TM automated homgeniser. The resulting mixture was then transferred to 50 mL polypropylene tubes and centrifuged (15 mins, 4000 g, 4 o C). The supernatant was then transferred to a clean polypropylene tube and dried under nitrogen with heat (55 o C) to an approximate volume of 1 mL. The extract is then transferred to 1.5 mL ultracentrifuge tubes and centrifuged at 14,000 g for 30 mins. The supernatant (200 µl) is then transferred directly to autosampler vials and 15 µl is injected onto the LC-MS/MS system.

Sample extraction and clean-up

Matrix-Extracted Calibration

Matrix extracted calibration curves were prepared and used for quantification. Control bovine muscle previously tested and found to contain no detectable residues was prepared as above. One muscle sample was used for each calibration standard level.

Samples were fortified at levels corresponding to 0, 1, 1.5, 2, 5 and 10 µg kg -1 by adding appropriate portions of the 100 ng mL -1 . After fortification, samples were held for 15 min prior to the extraction procedure as described above. Calibration curves were prepared by plotting the response factor (area of the analyte/area of internal

Method validation

The LC-MS/MS method was validated according to 2002/657/EC guidelines. LC-MS/MS identification criteria were verified throughout the validation study by monitoring relative retention times (see Table 1 for information with regard which internal standards are used with which analytes), signal to noise ratios (S/N) and ion ratios. Several method validation parameters were determined including linearity, specificity, recovery, precision (repeatability and within-laboratory reproducibility), decision limit CCα, and detection capability CCβ and measurement of uncertainty.

The decision limit (CCα) of the method was calculated according to the calibration curve procedure using the intercept (value of the response factor, y, where the concentration, x is equal to zero) and 2.33 times the standard error of the intercept for a set of data with 6 replicates at 3 levels. The detection capability (CCβ) was calculated by adding 1.64 times the standard error to the CCα. The specificity of the method was verified through the analysis of 20 blank muscle samples obtained from different sources. To test the linearity of the method, matrix extracted calibration curves were prepared on each of the three days to give 6-point calibration curves ranging from 0 to 10 µg kg -1 for each compound.

Since no certified reference materials were available for the analytes and matrix of interest, accuracy was determined from fortified negative samples as an alternative to trueness. The accuracy and precision were determined through the analysis of blank muscle samples fortified in 6 replicates at 1, 1.5, and 2 µg kg -1 on three separate days.

To evaluate the precision of the method, repeatability and within-laboratory The within-laboratory reproducibility results were obtained at the three validation levels on three separate days. The measurement of uncertainty was estimated by taking into account the within laboratory reproducibility over days 1, 2 and 3.

Results and Discussion

Preliminary experiments:-

The LC-MS/MS method was developed to provide confirmatory data for the analysis of bovine muscle for a wide number of steroidal type growth promoters. Electrospray ionisation (ESI) was chosen as the ionisation method as it provided the best responses on average for the wide array of compounds being investigated. Both positive and negative polarity ESI were tried for each of the analytes; some responded better in ESI-and others with ESI+, the final polarity mode chosen for each analyte can be seen in Table 1. The MS/MS fragmentation conditions were investigated and collision energies were optimised for each precursor-product ion transition, the fragmentation patterns observed were similar to those found in previously published work (Hewitt et al). For a method to be deemed confirmatory in accordance with Commission Decision 2002/657/EC four identification points are required. These identification points are obtained by monitoring one precursor ion and two transition products (corresponding to strong and weak ion transitions). The precursor and product transitions for each of the steroidal type compounds assayed for are listed in Table 1 as well as their associated collision energy. as has been seen before in numerous other publications. Extraction with organic solvent was then the initial step in this method. Various solvents (acetonitrile, methanol, diethylether, ethyl acetate) and combinations of these have been used

previously for the extraction of steroids from tissue samples. Due to the large number of steroids being examined in this study, a generic extraction solvent was required in order to minimize losses for some of the compounds. Acetonitrile was chosen initially for this purpose and initial studies indicated that the analytical recovery achieved was adequate to allow for the determination of these substances at and below1 µg kg -1 .

In relation to the purification step, again this needed to be quite generic in order to avoid analytical losses for some of the analytes. The possibility of using dispersive SPE with C18 sorbent material was investigated, however for some of the compounds examined; the analytical losses were too high to allow for determination at the level of interest. The addition of both sodium chloride and sodium sulphate during the initial extraction was also examined. Promising results were observed from the addition of sodium sulphate so an extraction procedure involving the addition of acetonitrile and sodium sulphate was chosen. A high speed centrifugation step was also implemented prior to LC-MS/MS analysis in order to aid the removal of particulates. The method proved to be quite robust with no deterioration in method performance being observed for all validation batches and survey samples.

The compounds were chromatographed on an Agilent Eclipse C 18 , 50 x 3mm, 1.8 µm column, the gradient profile was developed to permit rapid analysis but also to allow for the separation of the various isomers being examined which could not be separated by their m/z ratio; so needed to be separated chromatographically. Deuterated analogues of a number of the analytes were used as internal standards.

Validation study

Validation of the method was according to procedures described in Commission Decision 2002/657/EC covering specificity, calibration curve linearity, recovery (accuracy), precision, decision limit (CCα) and detection capability (CCβ).

Liquid Chromatography Tandem Mass Spectrometry

The criteria for relative retention times and ion ratios were examined for all samples and standards used for the validation study. The values for these were in agreement with the EU requirements for all the analytes investigated in this study. In terms of relative retention time, the analyte peaks in samples were found to be within 2.5% tolerance when compared with the matrix extracted standards. Furthermore two transition ions were monitored for each of the analytes. The most intense ion was used for quantitation. All ion ratios of samples were within the required tolerances as required by EU criteria when compared with standards used during the validation study.

F

Specificity

The technique of LC-MS/MS itself offers a high degree of selectivity and specificity.

To establish the selectivity/specificity of the method, a variety of muscle samples from different sources were fortified with the analytes of interest and their non-fortified equivalent samples were also analysed. No interfering peaks were observed at the retention time for any of the analytes.

Linearity of the response

The linearity of the chromatographic response was tested with matrix extracted curves using 6 calibration points in the concentration range of 0 to 10 µg kg -1 for each analyte. The regression coefficients (r 2 ) for all the calibration curves used in this study were ≥ ≥ ≥ ≥ 0.980.

Recovery

The recovery of the method was determined using muscle samples fortified at 1, 1.5 and 2 µg kg -1 . The resulting values were calculated from matrix extracted calibration curves. The mean corrected recovery (n = 18) for each of the analytes, determined over three separate assays is shown in Table 3 and ranged between 94.6 and 102.9 %.

Precision

The usefulness of suitable isotopically labelled internal standards is demonstrated in the excellent repeatability and within-laboratory reproducibility obtained for the compounds which have a corresponding deuterated analogue. No deuterated analogues were available for a number of these compounds and the precision values estimated for these were in the majority of cases worse than those observed for analytes with 13 deuterated analogues. The precision of the method, expressed as RSD values for the within-lab reproducibility at the three levels of fortification can be seen in Table 2.

CCα and CCβ

The decision limit is defined as the limit above which it can be concluded with an error probability of α, that a sample contains the analyte. In general, for prohibited substances an α value equal to 1 % is applied. The detection capability is the smallest content of the substance that may be detected, identified and quantified in a sample, with a statistical certainty of 1-β, were β = 5 %. CCα and CCβ were calculated using the intercept (value of the signal, y, were the concentration, x is equal to zero) and the standard error of the intercept for a set of data with 6 replicates at 3 levels on three separate days. Blank muscle was fortified at 1, 1.5 and 2 µg kg -1 . CCα is the concentration corresponding to the intercept + 2.33 times the standard error of the intercept. CCβ is the concentration corresponding to the signal at CCα + 1.64 times the standard error of the intercept (i.e. the intercept + 3.97 times that standard error of the intercept). The values determined for both CCα and CCβ can be seen in Table 2.

Measurement Uncertainty

The measurement of uncertainty was estimated by taking into account the within laboratory reproducibility over days 1, 2 and 3. The within laboratory reproducibility was multiplied by a coverage factor of three to give an overall figure for the uncertainty of the measurement. This approach of using the within laboratory reproducibility as a good estimator of measurement of uncertainty is taken from the SANCO/2004/2726rev1 document (SANCO 2004(SANCO /2726rev1 2004)). It recommends using the within laboratory reproducibility and using a coverage factor of 2.33 to estimate expanded uncertainty, however as it was felt that not all the environmental and other factors that could be varied over the course of the validation were examined, it was felt that a coverage factor of 2.33 may underestimate the true uncertainty of the method. So a value of 3 was chosen instead to give a more realistic value for the true uncertainty.

Surveillance of synthetic and natural growth promoters in imported bovine muscle samples

Criteria for confirmation of residue identity

The analytical method was designed to be used for both the screening and confirmation for 24 different natural and synthetic growth promoters. Several criteria based on those laid out by the European Union (2002/657/EC), all of which must be simultaneously satisfied were established for the confirmation of identity of a suspected residue. The first criterion is that both precursor-product ion transitions corresponding to the quantification and confirmation MRM transitions must contain peaks with a relative retention time (RRT) within 2.5% of the average RRT for the calibration standards under identical instrumental conditions. The second criterion is that the relative ratio of the peak areas in the precursor-product ion transitions corresponding to the quantification and confirmation MRM transitions must fall within an acceptable range of the average ratio observed for the calibration standards analysed with that batch of samples, the ranges for ion ratio deemed acceptable in accordance with 2002/657/EC are shown in Table 3. The calculated concentration of any analyte found in the survey samples must also be above the CCα determined for that compound for positive identification also.

Survey Results

A total of 102 bovine muscle samples imported from third countries into Ireland were collected over a period of three years (2007)(2008)(2009)(2010). No information was available on the breed, sex, age or species of the animals from which the muscle came from. Levels of progesterone were found in a number of samples at concentrations ranging between 0.28 µg kg -1 and 30.30 µg kg -1 . Levels of alpha and beta testosterone were also found in a number of samples at concentrations ranging between 0.22 to 8.63 µg kg -1 .and 0.16 to 2.08 µg kg -1 respectively. Table 4 presents a breakdown of the number of samples found to contain these three substances and also presents the concentration ranges found. There were a number of other samples where the criteria for positive identification were met for progesterone, alpha testosterone and beta testosterone, however the concentrations observed were less than the determined CCα values in each case. There were no observed levels above CCα for any of the synthetic steroidal growth promoters investigated as part of this study. However in one sample; the criteria for confirmation were met with respect to dexamethasone but the calculated concentration was well below the CCα value determined for this compound. Figure 3 shows chromatograms of one of the survey samples which was found to contain alpha testosterone at a concentration of 0.22 µg kg -1 and beta testosterone at a concentration of 0.37 µg kg -1 . The presence of progesterone, alpha testosterone and beta testosterone in tissue samples does not necessarily imply misuse of these substances, however one sample for contained a level of alpha testosterone approximately four times higher than the levels observed for the other 101 samples. The levels of progesterone found in muscle are in line with previous studies carried out in bovines with levels as high as 43.7 µg kg -1 being reported previously in bovine females [START_REF] Hartwig | Physiological quantities of naturally occurring steroid hormones (androgens and progestogens), their precursors and metabohtes in beef of differing sexual origin[END_REF] 1989). In relation to alpha testosterone; the highest level found in this study was 8.63 µg kg -1 , this is much higher than previously reported levels which ranged between 0.2-0.6 µg kg -1 [START_REF] Hartwig | Physiological quantities of naturally occurring steroid hormones (androgens and progestogens), their precursors and metabohtes in beef of differing sexual origin[END_REF].

Conclusions

The aim of this work was to develop and validate a rapid multi-residue LC-MS/MS confirmatory method that can simultaneously identify and quantify and quantify a large number of natural and prohibited steroidal type growth promoters in bovine muscle samples. This goal was successfully achieved with satisfactory results in terms of sensitivity, accuracy, linearity and repeatability being observed for all of the 24 analytes under investigation. The calculated CCα values were less than 0.28 µg kg -1 in all cases except for hexestrol where the calculated CCα value was 0.46 µg kg -1 .

The developed method was successfully used to carry out a small surveillance study into the presence of these 24 analytes in 102 samples of bovine muscle imported into Ireland from third countries over a period of three years. Application of the method demonstrated that it was capable of detecting the presence at very low levels of the naturally occurring steroids; progesterone, alpha testosterone and beta testosterone in a variety of these samples.

F

  C18 and Oasis MAX (Chursch et al 2008), Oasis HLB and aminopropyl (Schmidt et al 2009), graphitised carbon black and et al. 2009) C18 (Xu et al. 2006) and Strata X (Malone et al 2009).

  a function of analyte concentration (0 to 10 µg kg -1 ) and these were used to quantify samples.

  determined. Repeatability results were obtained at the three validation levels under identical conditions by the same operator on a particular day.

  out in this laboratory for the analysis of steroids in muscle involved a deconjugation step using Helix Pomatia. However the necessity of this step has been called into question for the analysis of tissue samples. Work carried out by others[START_REF] Hartmann | Simultaneous determination of anabolic and catabolic steroid hormones in meat by gas chromatography-mass spectrometry[END_REF][START_REF] Marchand | Ultra trace detection of a wide range of anabolic steroids in meat by gas chromatography coupled to mass spectrometry[END_REF], Yang et al. 2009) has shown that only a small percentage of the steroid exists as the glucuronide or sulphate form in muscle. It was decided to omit this step as part of this study so only the unconjugated steroid is being measured. Due to the absence of this deconjugation step, there was no need for the use of acetate buffers to bring the pH of the samples to approximately 5.2

  Number of survey samples which contained concentrations of Progesterone and alpha and beta Testosterone in the ranges 0.2-1 µg kg -1 , 1-5 µg kg -1 , 5-10 µg kg -1 , 10-20 µg kg -1 and >20 µg kg -1 .

  et al 1992). With regard beta testosterone; the highest concentration recorded in this study was 2.08 µg kg -1 unfortunately no data is available on the sex of the animals in this study but levels of beta testosterone in the range of 0.2-2.8 µg kg -1 have been

	observed previously for intact males (Hartwig et al 1997) All natural levels measured
	were below the tolerance value of 3.2 µg kg-1 for free testosterone in the muscle of
	fully grown cattle set by the Scientific Veterinary Committee (Scientific Veterinary
	F r o Committee
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	, Tsujioka

Table 1 :

 1 Figure 1 LC-MS/MS chromatograms (weak transition) of a blank muscle sample fortified with each of the analytes which ionize in ESI -at a level of 1 µg kg -1 . MS/MS Parameters for all investigated analytes.

	o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

* Indicates the most abundant transition; which is used for quantitation

Table 2 :

 2 CCα, CCβ, Inter-assay RSD (%), Measurement of Uncertainty (%) and Accuracy (%) values for each of the analytes investigated.

	Compound	CCα	CCβ	R.S.D.	M.U.	Accuracy
		(µ µ µ µg kg -1 )	(µ µ µ µg kg -1 )	(%)	(%)	(%)
	Medroxyprogesterone Acetate	0.12	0.20	4.4	13.2	101.7
	Melengestrol Acetate	0.13	0.22	4.8	14.4	100.1
	Megestrol Acetate	0.13	0.22	5.3	15.9	100.4
	Delmadinone Acetate	0.17	0.29	6.5	19.5	98.3
	Chlormadinone Acetate	0.17	0.30	6.4	19.2	99.5
	Methylboldenone	0.16	0.28	6.3	18.9	100.7
	Methyltestosterone	0.11	0.19	4.0	12	98.7
	Fluoxymesterone	0.20	0.34	7.4	22.2	101.7
	alpha Trenbolone	0.16	0.27	6.1	18.3	100.1
	beta Trenbolone	0.13	0.22	4.7	14.1	101.5
	16 beta hydroxy Stanozolol	0.17	0.29	7.3	21.9	99.3
	Progesterone	0.23	0.39	8.8	26.4	96.8
	alpha Testosterone	0.20	0.34	7.6	22.8	102.6
	beta Testosterone	0.16	0.27	5.7	17.1	99.7
	alpha Boldenone	0.15	0.25	6.6	19.8	100.6
	beta Boldenone	0.11	0.19	4.5	13.5	100.0
	alpha Nortestosterone	0.17	0.29	6.0	18.0	100.7
	beta Nortestosterone	0.16	0.28	6.2	18.6	102.2
	Taleranol	0.12	0.21	5.8	17.4	99.8
	Hexestrol	0.46	0.79	11.8	35.4	94.6
	Diethylstilbestrol	0.28	0.48	9.0	27.0	102.9
	Flumethasone	0.17	0.29	6.0	18	101.8
	Dexamethasone	0.13	0.21	4.8	14.4	100.1
	Zeranol	0.14	0.23	5.5	16.5	99.1

Table 3 :

 3 Maximum permitted tolerances for ion ratio values for LC-MS/MS analysis according to Commission Decision 2002/657/EC.

	Page 25 of 29

http://mc.manuscriptcentral.com/tfac Email: fac@tandf.co.

uk Food Additives and Contaminants

http://mc.manuscriptcentral.com/tfac Email: fac@tandf.co.ukFood Additives and Contaminants

Acknowledgements

This work was funded by The State Laboratory, Co. Kildare in Ireland. The authors would like to thank the staff at The State Laboratory, Co. Kildare for their practical assistance.


Ana L Cordeiro
email: cordeiro@ipfdd.de

Catharina Hippius
email: catharina.hippius@basf.com

Carsten Werner
email: werner@ipfdd.de

Immobilized enzymes affect biofilm formation

Keywords: bioactive surfaces, biofilm, cellulase, Pseudomonas aeruginosa, Staphylococcus epidermidis, subtilisin

The effect of the activity of immobilized enzymes on the initial attachment of pathogenic bacteria commonly associated with nosocomial infections (Pseudomonas aeruginosa and Staphylococcus epidermidis) was investigated. The proteolytic enzymes, subtilisin A and the glycoside hydrolase cellulose, were covalently attached onto poly(ethylene-alt-maleic) anhydride copolymer films. A comparison between active and heatinactivated surfaces showed that while the activity of immobilized cellulase reduced the attachment of S.

epidermidis by 67%, it had no effect on the attachment of P. aeruginosa. Immobilized subtilisin A had opposite effects: the active enzyme had no effect on the attachment of S. epidermidis but reduced the attachment of P. aeruginosa by 44%. The results suggest that different biomolecules are involved in the initial steps of attachment of different bacteria, and that the development of broad-spectrum antifouling enzymatic coatings will need to involve the co-immobilization of enzymes.

Introduction

In general, the exposure of a surface to aqueous environments results on its colonization by organisms. This phenomenon, known as biofouling, has severe negative impacts in industry and in medicine, for example by clogging membranes and filters, or causing infections. The use of enzymes to control biofouling has been receiving increasing attention [START_REF] Olsen | Enzyme-based antifouling coatings: a review[END_REF][START_REF] Kaplan | Therapeutic potential of biofilm-dispersing enzymes[END_REF][START_REF] Kristensen | Antifouling enzymes and the biochemistry of marine settlement[END_REF]Cordeiro and Werner in press 2011). Enzymes are currently available at affordable prices and are therefore viable for industrial use, and since they are biodegradable, they are also environmentally friendly. The most explored antifouling enzymatic modes of action include the degradation of cell wall components leading to cell lysis, and the degradation of extracellular polymer substances anchoring cells to surfaces (Cordeiro and Werner 2011).

Various laboratory-screening assays, studying the effects of enzyme formulations on the attachment of diverse foulers, identified the broad-spectrum serine protease subtilisin (the active component of Savinase, Alcalase and Esperase) as the most promising commercially available enzyme for antifouling strategies based on the degradation of adhesive biomolecules [START_REF] Pettitt | Activity of commercial enzymes on settlement and adhesion of cypris larvae of the barnacle Balanus amphitrite, spores of the green alga Ulva linza, and the diatom Navicula perminuta[END_REF][START_REF] Dobretsov | Inhibition of biofouling by marine microorganisms and their metabolites[END_REF]Leroy et al. 2008;Cordeiro and Werner 2011). Among the variety of enzymes tested to disperse biofilms, the glycoside hydrolase cellulase emerges as a common component of enzymatic solutions for removal and prevention of bacteria biofilm in industry. Studies showed that, although free cellulase affects Pseudomonas aeruginosa biofilm formation and mature biofilms [START_REF] Loiselle | The use of cellulase in inhibiting biofilm formation from organisms commonly found on medical implants[END_REF], crosslinked cellulase only partially inhibited P. aeruginosa biofilm formation [START_REF] Anderson | Use of enzymes to inhibit biofouling Proceedings 26th Annual Meeting of the Adhesion Society[END_REF].

Recently, we have developed a well-defined model system to investigate the influence of immobilized subtilisin A on the adhesion of major marine foulers (Tasso et al. 2009a;Tasso et al. 2009b). Our studies showed that the adhesion strength of Ulva linza spores and Navicula perminuta cells decreased in the presence of the active enzyme, and that the effect of the enzyme was significantly enhanced when the enzyme was surfaceconfined (Tasso et al. 2009b), revealing the important role of enzyme location at the surface-fouler interface on antifouling efficacy.

In this paper, we evaluate the impact of immobilized biofilm-degrading enzymes on the initial attachment of P. aeruginosa and Staphylococcus epidermidis, bacteria commonly associated with hospital-acquired infections.

To achieve this purpose, subtilisin A and cellulase were covalently immobilized onto poly(ethylene-alt-maleic) anhydride copolymer films.

Materials and methods

Preparation of bioactive layers

The maleic anhydride (MA) copolymer films were prepared as described elsewhere [START_REF] Pompe | Maleic anhydride copolymers -A versatile platform for molecular biosurface engineering[END_REF]. Briefly, silicon wafers or glass slides were freshly oxidized in a mixture of aqueous solutions of NH 3 and H 2 O 2 and thereafter surface-modified with 3-aminopropyltriethoxysilane (APTES) (Fluka, Germany) by immersing the freshly oxidized wafers in a solution of APTES in a 9:1 (v/v) mixture of 2-propanol and water followed by extensive rinsing with 2-propanol. Stable covalent binding was achieved by annealing at 120 °C for 1 h.

Copolymer solutions of poly(octadecene-alt-maleic anhydride) (POMA) MW 30 000 -50 000 (Polysciences Inc., USA) and poly(ethylene-alt-maleic anhydride) (PEMA) MW 125 000 (Aldrich, Germany) were prepared by dissolving POMA in THF at 0.08% (w/v), and PEMA at 0.15% (w/v) in THF/acetone (2:1 v/v). The copolymer solutions were spin-coated (Suess Microtec, Germany) onto the amino-functionalized wafers. Stable covalent binding of the polymer films was achieved by annealing at 120 °C to generate imide bonds with the aminosilane [START_REF] Pompe | Maleic anhydride copolymers -A versatile platform for molecular biosurface engineering[END_REF]. The annealed samples were extensively rinsed to wash out polymer excess.

The properties of the films were investigated using a wide range of techniques, details can be found elsewhere [START_REF] Pompe | Maleic anhydride copolymers -A versatile platform for molecular biosurface engineering[END_REF][START_REF] Pompe | Functional films of maleic anhydride copolymers under physiological conditions[END_REF][START_REF] Osaki | Ionization characteristics and structural transitions of alternating maleic acid copolymer films[END_REF].

Reactive immobilization of enzymes onto MA copolymer films was achieved by exposing the copolymer layers to enzyme solutions of variable concentration, as described in (Cordeiro et al. 2011;Tasso et al. 2009a).

Immediately prior to enzyme immobilization, the anhydride moieties of the copolymers were regenerated by annealing at 120 °C for 2 h. The cellulase solutions were prepared by dissolving cellulase in 50 mM sodium citrate buffer pH 4.0. Subtilisin A solutions were prepared by dissolving the enzyme in PBS modified by the addition of CaCl 2 and NaCl at 2 mM and 0.1 M, respectively; the pH of the solution was adjusted to 8.6 with NaOH (Tasso et al. 2009a). The immobilization was performed using in-house constructed immobilization chambers. After overnight exposure of the copolymer films to the enzyme solution the samples were rinsed 10 times with deionized water to remove unbound protein (Tasso et al. 2009a).

Characterization of bioactive layers

The immobilization and characterization of subtilisin A onto MA copolymer films are described in detail elsewhere (Tasso et al. 2009a;2009b).

The dry thickness of the cellulase layers prepared onto MA copolymer films were determined by singlewavelength ellipsometry (Sentech, Germany) at an angle of incidence of 70° at 632.8 nm. Upon determination of the thickness of all layers below the protein layer, a five-model approximation (Si/SiO 2 /silane + MA film/protein/air) was used to determine the thickness of the protein layer. The refractive index for the MA copolymer film was taken from [START_REF] Osaki | Ionization characteristics and structural transitions of alternating maleic acid copolymer films[END_REF]. The refractive index of the protein layers was maintained fixed for all protein-containing surfaces at 1.465 [START_REF] Fischer | The ability of surface characteristics of materials to trigger leukocyte tissue factor expression[END_REF]. All ellipsometry measurements were carried out at 20 °C and 55% air humidity.

The activity of immobilized cellulase was determined through the detection of released reducing sugars from carboxymethylcellulose (CMC) using 3,5-dinitrosalicyclic acid (DNS) adapting the procedures described in [START_REF] Ghose | Measurement of cellulase activities[END_REF]. Briefly, substrate solutions consisting of 2% w/v CMC in 0.05 M sodium citrate buffer pH 4.8 were prepared by stirring the solution overnight at room temperature and subsequently stored at 4 °C. Activity assays were performed by pre-incubation of the CMC substrate at 50 °C and by immersion of the immobilization chambers in a water bath kept at 50 °C. Upon temperature equilibration, the substrate solution was added at constant volume to each surface. At intervals, 50 l was extracted and added to respective wells in PCR plates already containing 150 l DNS reagent. (DNS reagent was prepared by dissolving 0.5 g 3,5-dinitrosalicyclic acid in 20 ml 1 M NaOH and 15 ml distilled water after which 15 g potassium sodium tartrate was added and the solution was made up to 50 ml with distilled water.) Using a thermal cycler, the PCR plates were immediately heated up to 99 °C, kept at this temperature for 5 min and then cooled down to room temperature. Aliquots were extracted, transferred to microtiter plates and the volume made up to 200 l with distilled water. Absorbance was measured at 540 nm. Blanks consisting of the substrate solutions not contacting active enzyme were used. A calibration curve using glucose was prepared to allow for the conversion of the measured absorbance to the equivalent concentration of released glucose.

Experiments with bacteria

The antifouling potential of the activity of immobilized enzymes was investigated using biofilms of single bacterial species (Pseudomonas aeruginosa (DSMZ 939; ATCC 15442) and Staphylococcus epidermidis (DSMZ 1798; ATCC 12228)). The bacterial were kept on agar plates (TSBY and LB agar for S. epidermidis and P. aeruginosa, respectively) at 2 °C in the dark.

Test samples consisting of PEMA copolymer films decorated with active enzyme were compared with control samples consisting of PEMA copolymer films exposed overnight to enzyme buffer (i.e. polymer coatings not containing enzyme), and with control samples consisting of PEMA copolymer films containing heatdenatured enzyme. The heat-denatured subtilisin A controls were prepared by heating immobilized subtilisin A layers to 120 °C for 45 min as described in (Tasso et al. 2009b). Inactivated cellulase layers were prepared by autoclaving correspondent active layers for 20 min at 120 °C. Immediately after enzyme immobilization, test and control samples (in triplicate) were incubated with bacteria suspensions of approximately 10 6 c.f.u./ml at 37 °C for 24 h (humidity >90%) under static conditions using sterile in-house constructed chambers (similar to those used for enzyme immobilization). Suspensions of S. epidermidis and P. aeruginosa were prepared in 5% TSBY medium for S. epidermidis, and in 10% LB medium for P. aeruginosa. Biofilm growth conditions were optimized with respect to bacterial cell and media concentrations, as well as incubation times for each strain, to ensure reproducible and homogenous biofilm formation (data not shown). After incubation, the slides were removed from the chambers and washed three times to remove loosely bound bacteria by immersion in 50 ml saline solutions in Falcon conical centrifuge tubes. Samples were then transferred to tubes containing 50 ml saline and submitted to 1 min vortexing 5 min ultrasonication (VWR ultrasonication cleaner USC300D, power level 5) and 1 min vortexing. The cells recovered from the surfaces after washing and ultrasonication (i.e. sessile bacteria) were plated in serial dilutions for conventional viable cell count analysis. Biological data were tested towards significance by pairwise comparisons within the analysis applying Students t-Test.

Results and discussion

Bioactive layers

The strategy used to graft PEMA and POMA copolymer layers onto SiO 2 surfaces, and to covalently immobilize the enzymes, is outlined in Figure 1. For the covalent binding of MA copolymers, primary amino functionalities were introduced on the glass surfaces by binding an APTES precursor layer. Upon elevation of temperature, the introduced primary amines reacted with the anhydride group of the MA copolymers forming imide bonds.

Subtilisin A and cellulase were then covalently immobilized onto PEMA and POMA copolymer films through the spontaneous reaction of free, non-reacted anhydride groups with primary amines of the enzymes. The variation of the concentration of the immobilization solutions further allowed obtaining enzyme-decorated films of distinct properties.

INSERT FIGURE 1

MA copolymer films covered with subtilisin A, were extensively characterized concerning immobilized amount (determined by ellipsometry, amino acid analysis using HPLC and confocal laser scanning microscopy using fluorescently labeled enzyme) and activity (determined by following the degradation of Suc-AAPF-pNA by absorbance spectroscopy at 405 nm). Further details can be found in Tasso et al. (2009a). Briefly, the amount and activity of immobilized subtilisin A was dependent on the concentration of the immobilization solution and of the copolymer film. Highly swelling PEMA copolymer layers enabled the immobilization of up to five times more enzyme than hydrophobic POMA copolymer films.

Since cellulase consists of a complex mixture of more than three protein components (for details see [START_REF] Beldman | The cellulase of Trichoderma viride[END_REF], the quantification of the immobilized amount was not possible by a.a. analysis using HPLC. Therefore, the amount of immobilized cellulase was estimated by measuring the thickness of the immobilized protein layer by ellipsometry (Fig. 2). The results show a thinner enzyme layer on POMA copolymer films as compared with PEMA. The compact, non-swelling properties of the POMA copolymer films [START_REF] Pompe | Functional films of maleic anhydride copolymers under physiological conditions[END_REF]) possibly triggers conformational changes and spreading of the proteins onto the surface as a result of hydrophobic interactions between the surface and the hydrophobic moieties of the proteins.

Immobilization onto PEMA copolymer films enabled the generation of enzyme layers up to 10 times thicker than onto POMA layers, with minor dependence on the concentration of the immobilization solution (within the range tested). The hydrophilic, highly swelling, 3-D-like PEMA copolymer film possibly enables the penetration of the enzyme into the swollen layer resulting in higher loading. In addition, cellulose-containing PEMA copolymer films prepared with enzyme concentrations below 2 mg/ml were slightly inhomogeneous (results not shown). The accelerated transport of the proteins to the surface resultant from more concentrated immobilization solutions possibly favors cellulase immobilization onto hydrophilic polymer film.

INSERT FIGURE 2

The activity of the cellulase layers was determined by following the degradation of CMC at 50 °C using DNS to detect the amount of released reducing sugars [START_REF] Ghose | Measurement of cellulase activities[END_REF]. The initial activity, expressed as the concentration of released glucose equivalent per unit time, was determined through the slope of the linear region of the reaction progress curve. The initial activity of the cellulase layers in dependence of the copolymer film and of the concentration of the enzyme immobilization solution is presented in Fig. 3.

INSERT FIGURE 3

The results show that for immobilization of cellulase below 10 mg/ml, the activity of the immobilized enzyme layer for the different MA copolymer films was similar. Since thicker enzyme layers were measured onto PEMA copolymer films as compared with POMA copolymer films in this concentration range, the results suggest a more efficient immobilization onto the hydrophobic coating (i.e. a higher amount of active enzyme per total immobilized amount). The low immobilization efficiency onto PEMA copolymer films may be attributed to limited substrate diffusion and/or to restrictions in structural transitions of the enzyme within the multilayer-like packed enzyme, as compared with POMA copolymer films. For immobilization solutions higher than 10 mg/ml, the immobilization of cellulase onto PEMA copolymer films enables generating cellulase films with approx. double the activity than onto POMA copolymer films. The higher concentration of the immobilization solution results in an accelerated transport of the enzyme to the polymer surface, possibly resulting in a conservation of enzyme conformation and activity.

Bacteria attachment

The effect of the activity of immobilized cellulase and of immobilized subtilisin A on the initial attachment of two bacteria species associated with hospital-acquired infections was investigated by comparing active enzyme containing PEMA copolymer coatings of highest activity (see Table 1 for coating properties) with respective heat-inactivated controls, and with bare PEMA copolymer films.

INSERT TABLE 1

The effects of immobilized cellulase and subtilisin A on the attachment of Staphylococcus epidermidis are presented on Fig. 4. While cell attachment was significantly reduced onto active cellulase layers, with 67% less cells attaching to the active surface as compared with respective denatured control, no significant differences were observed between active and inactive subtilisin A layers. The results suggest that proteins are playing a minor role on the initial attachment of S. epidermidis.

INSERT FIGURE 4

The effects of the activity of immobilized cellulase and of immobilized subtilisin A on the initial attachment of the pathogen P. aeruginosa is presented on Fig. 5. The results show that although immobilized subtilisin A reduced the attachment of P. aeruginosa by 44% when compared with its denatured control, immobilized cellulase had no significant effect on the attachment of P. aeruginosa when comparing active and denatured samples. In a first analysis, the results do not agree with the antifouling effect of free and immobilized cellulase towards P. aeruginosa as reported previously [START_REF] Anderson | Use of enzymes to inhibit biofouling Proceedings 26th Annual Meeting of the Adhesion Society[END_REF][START_REF] Loiselle | The use of cellulase in inhibiting biofilm formation from organisms commonly found on medical implants[END_REF]. The different effects observed might be explained by the different experimental conditions used, as for example surface enzymatic activity, the different bacteria strains and biofilm growth conditions, as well as the different methodologies used for biofilm quantification.

INSERT FIGURE 5

The results suggest that different biomolecules are involved in the initial steps of attachment of the different bacteria species tested. Proteins may be playing a larger role in the initial attachment of P. aeruginosa to surfaces as compared with S. epidermidis, for which polysaccharides seem to play the major role.

The mechanisms by which enzymes may exert an antifouling effect include the degradation of molecules anchoring the cells to the surfaces (e.g. adhesins), and the degradation of the extracellular matrix thereby weakening the biofilm structure and leading to dispersion or to the facilitated penetration of other antimicrobial compounds into the biofilm (Cordeiro and Werner, in press 2011). Previous studies showed that free subtilisin was more effective at inhibiting adhesion than at enabling detachment of the marine bacteria Pseudoalteromonas sp. D41 (Leroy et al. 2008b). In contrast, the enzyme inhibitory effects were comparable to the detachment effects for a multispecies biofilm [START_REF] Hangler | The serine protease Esperase HPF inhibits the formation of multispecies biofilm[END_REF] suggesting that, in this case, the enzyme effectively degrades both protein-based adhesives and proteins contained in the matrix [START_REF] Hangler | The serine protease Esperase HPF inhibits the formation of multispecies biofilm[END_REF]. Recent work also showed that differences in chemical composition of the EPS are reflected on the vulnerability of biofilms to enzymatic treatment [START_REF] Lequette | Using enzymes to remove biofilms of bacterial isolates sampled in the food-industry[END_REF][START_REF] Chaignon | Susceptibility of staphylococcal biofilms to enzymatic treatments depends on their chemical composition[END_REF][START_REF] Augustin | Assessment of enzymatic cleaning agents and disinfectants against bacterial biofilms[END_REF]. For example, while staphylococcal biofilms producing substantial amounts of the polysaccharide adhesin poly-N-acetylglucosamine (PNAG) could be degraded by Dispersin B, biofilms not producing PNAG could not be dispersed by Dispersin B but by proteases [START_REF] Chaignon | Susceptibility of staphylococcal biofilms to enzymatic treatments depends on their chemical composition[END_REF].

Our results suggest that proteinaceous adhesives play a larger role in the attachment of P. aeruginosa than in the attachment of S. epidermidis. Further studies on the effects of immobilized cellulase and subtilisin A on the adhesion strength of the attached bacteria (i.e. testing the enzyme dispersal properties), and on possible synergetic effects of the activity of the immobilized enzyme and other antimicrobial agents (e.g. antibiotics) would be particularly important when considering the design of novel antimicrobial bioactive surfaces.

Altogether the results demonstrate that due to the complexity and variability of microbial attachment mechanisms and of EPS composition, the production of enzyme-based antimicrobial coatings targeting the degradation of adhesive components should involve the co-immobilization of different enzymes to achieve high and broad-spectrum antifouling efficacy. Table 1 Description

Tables

Figure captions Fig. 1

 captions1 Figure captions

Fig. 2

 2 Fig. 2 Thickness of the immobilized cellulase layer onto PEMA () and POMA () copolymer films in dependence of the concentration of the immobilization solution determined by ellipsometry. Enzyme immobilization performed overnight at RT in 50 mM sodium citrate pH 4.0. Error bars = ± SD.

Fig. 3

 3 Fig. 3 Initial activity (expressed as the amount of released glucose equivalent per min) of immobilized cellulase onto PEMA () and POMA () copolymer films in dependence of the concentration of cellulase used for immobilization. Enzyme immobilization performed overnight at RT in 50 mM sodium citrate pH 4.0. Activity measured by detection of released sugars (using DNS) from 2% (wt/v) CMC in 50 mM sodium citrate buffer pH 4.8 at 50 °C. Error bars = ± SD.

Fig. 4

 4 Fig. 4 Surface adherent Staphylococcus epidermidis cells onto PEMA copolymer films containing no, active, or heat-denatured (a) cellulase and (b) subtilisin A after 24 h incubation at 37 °C. Error bars = ± SD.

Fig. 5

 5 Fig. 5 Surface adherent Pseudomonas aeruginosa cells onto PEMA copolymer films containing no, active or heat-denatured (a) cellulase, and (b) subtilisin A after 24 h incubation at 37 °C. Error bars = ± SD.

Acknowledgments

The authors are grateful to Tina Lenk and Nelly Rein (Leibniz Institute of Polymer Research Dresden, Germany), and to Irina Elert (BASF SE, Ludwigshafen, Germany) for technical support.


Zoltan Dombovari

David A W Barton
email: david.barton@bristol.ac.uk

R Eddie Wilson
email: re.wilson@bristol.ac.uk

Gabor Stepan
email: stepan@mm.bme.hu

On the Global Dynamics of Chatter in the Orthogonal Cutting Model

Keywords: orthogonal cutting, delay differential equation, differential algebraic equation, non-smooth, chatter, fly-over

The large-amplitude motions of a one degree-of-freedom model of orthogonal cutting are analysed. The model takes the form of a delay differential equation which is non-smooth at the instant at which the tool loses contact with the workpiece, and which is coupled to an algebraic equation that stores the profile of the cut surface whilst the tool is not in contact. This system is approximated by a smooth delay differential equation without algebraic effects which is analysed with numerical continuation software. The grazing bifurcation that defines the onset of chattering motion is thus analysed as are secondary (period-doubling etc.) bifurcations of chattering orbits, and convergence of the bifurcation diagrams is established in the vanishing limit of the smoothing parameters. The bifurcation diagrams of the smoothed system are then compared with initial value simulations of the full non-smooth delay differential algebraic equation. These simulations mostly validate the smoothing technique and show in detail how chaotic chattering dynamics emerge from the non-smooth bifurcations of periodic orbits.

Introduction

In the machining industry it is well known that at high material removal rates, turning, milling and drilling processes are subject to chattering motions in which the tool repeatedly loses and re-establishes contact with the workpiece. Chattering results in a very poor quality finish, see figure 1(a), so in practice the need to avoid it limits the material removal rate and places bounds on the technological parameters such as the spindle speed and the depth of cut (also known as the chip width). These bounds are in addition to those due to the power and torque characteristics of the machine.

The general aim of this paper is to understand in more detail the dynamical mechanisms involved in the onset of chattering in turning processes, where a quasi-stationary tool cuts into a rapidly rotating workpiece. To simplify matters, we restrict our analysis to the special case of orthogonal cutting, depicted in figure 2.

In turning processes, the cutting force between the tool and the workpiece is a function of the chip thickness, that is, the difference between the tool position and the surface position. The surface position was determined one (or more) revolutions earlier by the tool's past motion. Turning processes are thus regenerative [START_REF] Tlusty | Self-excited vibrations on machine tools[END_REF][START_REF] Tobias | Machine-tool Vibration[END_REF] in the sense that the tool's motion is forced by its own history, leading to models based on delay differential equations (DDEs). These models admit equilibria corresponding to steady cutting, where the chip thickness is independent of time, but the stability of this regime is not assured. By analysing when steady cutting is linearly unstable, we may thus identify combinations of technological parameters which the machinist should avoid [START_REF] Altintas | Analytical prediction of stability lobes in milling[END_REF][START_REF] Budak | Analytical prediction of chatter stability in milling-part i: General formulation[END_REF][START_REF] Merdol | Multi frequency solution of chatter stability for low immersion milling[END_REF][START_REF] Insperger | Semi-discretization method for delayed systems[END_REF][START_REF] Bayly | Stability of interrupted cutting by temporal finite element analysis[END_REF][START_REF] Zatarain | Analysis of the influence of mill helix angle on chatter stability[END_REF][START_REF] Insperger | Criticality of Hopf bifurcation in state-dependent delay model of turning processes[END_REF][START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF].

In practice, the stability of the process is complicated by nonlinear effects, and the typical situation is indicated by the one-dimensional bifurcation diagram in figure 1(b), where we fix the spindle speed and analyse the dynamics of the tool as the chip width is varied. Here x := max t x(t)min t x(t) denotes the peak-to-peak magnitude of the tool's motion x(t), so that on branch A, where x = 0, we have linearly stable steady cutting in which the chip thickness is constant. As the chip width is increased, steady operation becomes linearly unstable at the Hopf bifurcation point B. However, normal form computations [START_REF] Stépán | Nonlinear regenerative machine tool vibrations[END_REF][START_REF] Szalai | Global dynamics of low immersion highspeed milling[END_REF] and measurements [START_REF] Shi | Theory of finite-amplitude machine-tool instability[END_REF][START_REF] Kalmár-Nagy | Experimental and analytical investigation of the subcritical instability in metal cutting[END_REF] show that the bifurcation is subcritical, indeed robustly so over a wide range of parameters and cutting force characteristics [START_REF] Dombóvári | Estimates of the bistable region in metal cutting[END_REF]. Hence the branch C of periodic orbits that emanates 'bends back' and the periodic orbits on it are themselves unstable and hence are not observable in experiments.

However, we know from experiments [START_REF] Shi | Theory of finite-amplitude machine-tool instability[END_REF][START_REF] Stepan | Modelling nonlinear regenerative effects in metal cutting[END_REF] that complicated large amplitude chattering motions E (figure 1(b)) are possible at depths of cut less than that of the Hopf bifurcation point. Hence in regime 2, there is bistability in that sufficiently large perturbations to steady cutting can result in sustained chattering oscillations [START_REF] Shi | Theory of finite-amplitude machine-tool instability[END_REF]. In contrast, in regime 1, only steady cutting is possible and in regime 3, only chattering is possible. Qualitatively this type of behaviour is seen in many different engineering applications, for example, pipe-flow [START_REF] Willis | Experimental and theoretical progress in pipe flow transition[END_REF].

The focus of this paper is the vicinity of the point D marked in figure 1(b). Specifically, we shall examine how the unsta- ble branch C connects to stable chattering motions in region E, and how the complexity of these motions develops through secondary bifurcations. The simplest mechanism for a periodic branch to regain stability is via a cyclic fold / SNLS (Saddle Node of a Limit Cycle) bifurcation, and we show that indeed this is the mechanism at D, albeit a non-smooth variant where the branch C turns back at the instant its orbits attain sufficient amplitude for the tool to leave the workpiece. This nonsmooth fold point has a rather intricate structure of accumulating branches which may be detected by following through sequences of further (non-smooth) fold and (smooth) period doubling bifurcations. See the monograph [START_REF] Bernardo | Piecewisesmooth dynamical systems: theory and applications[END_REF] for a recent account of non-smooth bifurcation theory.

To analyse these solution structures, our chief weapon is numerical continuation software, which is a well established technique for investigating dynamical systems. The key idea is that one may trace out branches of solutions (which may be equilibria or periodic orbits) as a parameter is varied without resorting to repeated solution of the initial value problem. In particular, branches of unstable solutions may be followed which subsequently bifurcate to produce stable dynamics, so that competing stable behaviours can be identified efficiently without an exhaustive search of the initial value and parameter spaces. The theoretical foundations of numerical continuation were developed in the 1980s (see for example [START_REF] Doedel | Numerical analysis and control of bifurcation problems, part I[END_REF][START_REF] Doedel | Numerical analysis and control of bifurcation problems, part II[END_REF][START_REF] Seydel | Practical bifurcation and stability analysis[END_REF] for an overview), and well-known packages include AUTO [START_REF] Doedel | AUTO 97: continuation and bifurcation software for ordinary differential equations[END_REF] and MatCont [START_REF] Dhooge | MATCONT and CL_MATCONT: Continuation toolboxes in Matlab[END_REF] for ordinary differential equations (ODEs). The analysis of DDEs, such as arise from regenerative processes, is complicated by the infinite dimensional state space, but these problems may be tackled by more sophisticated packages such as DDE-BIFTOOL [START_REF] Engelborghs | DDE-BIFTOOL v. 2.00: a Matlab package for bifurcation analysis of delay differential equations[END_REF] (which we use) and PDDE-CONT [START_REF] Szalai | PDDE-CONT: a continuation and bifurcation software for delay-differential equations[END_REF].

Unfortunately, numerical continuation depends upon continuous differentiability of the model with respect to its solutions and parameters, and our model for cutting is non-smooth at the instant when the tool loses contact with the workpiece. Furthermore, the tool position and its history are not sufficient to close the initial value problem for chattering motions. This is because when the tool is flying (i.e., has lost contact with the workpiece) a separate variable must be maintained to store the surface height of the workpiece for future computation after the tool has landed. This switch gives rise to a combined differential-algebraic system.

The key achievements of this paper are thus two-fold: 1) a smooth approximation to the cutting model, with no algebraic effects, is developed which may then be analysed with the DDE-BIFTOOL package, and 2) the bifurcation structure of the smooth approximation is elucidated and related to that of the non-smooth model, and, in particular, we show that the onset of chattering is via a type of non-smooth fold point.

The remainder of the paper is structured as follows. In section 2 we develop and parametrise a regenerative one degreeof-freedom model for orthogonal cutting, which is generalised to incorporate non-smooth and algebraic effects due to the tool leaving the workpiece. In section 3 we then describe how the model may be smoothed for computational purposes, and we show convergence of the bifurcation diagrams in the limit of vanishing smoothing parameters that demonstrate unequivocally that point D is a non-smooth cyclic fold. Then in section 4, we probe the secondary bifurcation sequences and the stable chattering behaviour in more detail and show how these correspond to numerical simulations of the initial value problem. Then section 5 presents conclusions.

Model description

In this section, a simple model is developed for the orthogonal cutting process that can demonstrate chattering behaviour. The model is formulated as a delay differential algebraic equation (DDAE) and includes the regenerative effect of turning processes and the non-smoothness when contact between the cutting tool and the workpiece is lost.

Permanent contact model

We begin by modelling the simple case where the tool is in permanent contact with the workpiece [START_REF] Dombóvári | Estimates of the bistable region in metal cutting[END_REF]. As we have described, turning processes have a so-called regenerative effect or memory property: when the tool is in permanent contact with 2) where the forces on the cutting tool are proportional to the chip thickness h(t) = q(t -τ)q(t) + h 0 , where τ = 2π/Ω is the period of revolution.

q 0 q t ( ) F W q x k b workpiece machine chip t ) ( - t q ( ) b F m q
the workpiece, the thickness of the chip that is cut is the difference between the current position of the cutting tool and the position of the cutting tool one revolution earlier. Consequently, this introduces a delay of time τ = 2π/Ω into the system, where Ω is the spindle speed. If q(t) denotes the displacement of the tool orthogonal to the workpiece, then the instantaneous chip thickness h(t) is given by

h(t) = q(t -τ) -q(t) + h 0 = Δh(t) + h 0 , (1)
where h 0 is the desired steady-state chip thickness (proportional to the feed rate), see figure 2(b).

We consider a one degree-of-freedom model of the motion of the cutting tool-toolholder-machine system as illustrated in figure 2(b). In the direction perpendicular to the workpiece the equation of motion is

q(t) + 2κω n q(t) + ω 2 n q(t) = 1 m F q , (2)
where q(t) is the position of the cutting tool relative to the tool holder that provides the constant feed, which is just h 0 for one revolution of the workpiece. The system parameters are the natural angular velocity ω n , the damping ratio κ and the modal mass m. Here F q is the orthogonal component of the cutting force which must be determined empirically. In practice, we find F q is linearly proportional to the chip width (depth of cut) w, but depends nonlinearly on the chip thickness h(t). Via (1), it follows that ( 2) is a DDE. Throughout the paper we fix κ = 0.04 as a representative damping ratio. The classical approximation to F q is given by the threequarter rule F q ∼ wh 3/4 , see [START_REF] Kienzle | Spezifische schnittkräfte bei der metallbearbeitung[END_REF], which has infinite gradient at h = 0 when the tool leaves the workpiece. To avoid this degeneracy, we use instead the approximation due to Shi and Tobias [START_REF] Shi | Theory of finite-amplitude machine-tool instability[END_REF], which takes the cubic polynomial form

F q (h(t)) = w(ρ 1 h(t) + ρ 2 h 2 (t) + ρ 3 h 3 (t)), (3)
with typical measured parameter values ρ 1 = 6109.6 N mm -2 , ρ 2 = -54141.6 N mm -3 and ρ 3 = 203769 N mm -4 for steel. We consider tool motion x(t) relative to the equilibrium position q 0 = F q (h 0 )/mω 2 n , and rescale time and displacement us-ing the natural angular velocity ω n and the desired chip thickness h 0 respectively, so that

t = ω n t, x = x/h 0 , w = w ρ 1 + 2ρ 2 h 0 + 3ρ 3 h 2 0 mω 2 n , τ = ω n τ, Ω = 2π τ = Ω/ω n , η 2 = h 0 (ρ 2 + 3ρ 3 h 0 ) ρ 1 + 2ρ 2 h 0 + 3ρ 3 h 2 0 and η 3 = h 2 0 ρ 3 ρ 1 + 2ρ 2 h 0 + 3ρ 3 h 2 0 , (4)
where the dimensionless cutting force parameters become η 2 = -0.20015 and η 3 = 0.01946. (These parameters are typical for finishing-type operations where surface quality is important.) Thus (2) may be rewritten in the canonical non-dimensional form [START_REF] Merdol | Multi frequency solution of chatter stability for low immersion milling[END_REF] where the tildes are dropped for clarity. Here Δh(t

x (t) + 2κx (t) + x(t) = F(Δh(t)) = w(Δh(t) + η 2 Δh 2 (t) + η 3 Δh 3 (t)),
) := h(t) - 1 = x(t -τ) -x(t) is the (non-dimensional) instantaneous chip thickness.

Non-smooth intermittent contact model

We now consider how to augment model ( 5) to include motions where the tool loses contact with the workpiece. At the instant that contact is first lost, the chip thickness attains zero and the cutting force is 'switched off' until the tool regains contact. In consequence, the definition of the cutting force is augmented with F(h) = 0 for h < 0, and thus it is not differentiable at h = 0.

Furthermore, since the original equation of motion (5) tracks only the position of the cutting tool and not the surface of the workpiece, equation (1) defining the chip thickness becomes invalid once contact is lost. Thus we introduce the surface height function χ(t) defined by the recurrence relation

χ(t) = x(t), if cutting, χ(t -τ) + 1, if flying, (6)
where x(t) is the position of the cutting tool from [START_REF] Merdol | Multi frequency solution of chatter stability for low immersion milling[END_REF]. Clearly, when the cutting tool is in contact with the surface, the positions of the two are identical. However, when the cutting tool is not in contact, the surface height is simply the same as it was one revolution earlier, corrected by the (non-dimensionalised) feed-rate. Consequently, the new expression for chip thickness becomes

h(t) = χ(t -τ) -x(t) + 1 = Δh(t) + 1. (7)
In summary, the final non-smooth model is a DDAE defined by ( 5), ( 6) and [START_REF] Bayly | Stability of interrupted cutting by temporal finite element analysis[END_REF]. For convenience we define the cutting and flying vector-fields S c and S f respectively in the form

S c (t) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ x (t) F(Δh(t)) -2κx (t) -x(t) x(t) -χ(t) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (8)
and

S f (t) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ x (t) -2κx (t) -x(t) χ(t -τ) + 1 -χ(t) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (9)
so that the full model may be summarised by

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ x (t) x (t) 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = H(h(t))S c (t) + [1 -H(h(t))]S f (t), ( 10
)
where H is the Heaviside step function. The (infinite dimensional) phase-space of ( 10) is given by (x t (s), x t (s), χ t (s)) for s ∈ [-τ, 0], where (•) t (s) = (•)(t + s), and consequently solutions of the model must be projected onto a two-dimensional surface in order to be visualised. In the analysis that follows we shall speak of the switching surface, defined to be the set of all states where the switching condition h(t) = 0 is satisfied, and thus where [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF] switches between the cutting and flying vector fields.

Smoothed equations

We now approximate the non-smooth system (10) of DDAE with a smoothed delay differential equation DDE approximation, so that even chattering motions may be analysed by the standard numerical continuation package DDE-BIFTOOL. We then demonstrate continuation computations and study their convergence in the smoothing parameters. In particular, we study the nature of the onset of chattering motions (point D in figure 1(b)) and verify that this a type of non-smooth fold.

Smoothed intermittent contact model

Firstly, in [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF], we may approximate the Heaviside step function with the smooth function

H(h(t)) ≈ H ε1 (h(t)) := 1 2 1 + tanh h(t) ε 1 , (11)
where ε 1 > 0 is a smoothing parameter, so that we expect to recover the full non-smooth dynamics in the limit ε 1 → 0. Secondly, we must address the more serious difficulty in [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF], namely the algebraic recursion, which cannot be tackled by standard numerical continuation packages. Our approach is to approximate (6) by a differential equation by introducing a further smoothing in the form of the singular perturbation ε 2 χ (t) (ε 2 > 0) in the third component of [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF]. The fully smoothed DDE approximation thus takes the form

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ x (t) x (t) ε 2 χ (t) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = H ε1 (h(t))S c (t) + [1 -H ε1 (h(t))]S f (t), ( 12
)
where χ is now a fast variable that tracks either x(t) or χ(t-τ)+1 with the characteristic time scale ε 2 .

The choice of the smoothing parameters ε 1 and ε 2 is critical and we experiment with them later in section 3.3. In particular, they must be chosen sufficiently small so as to avoid introducing dynamical artefacts which are solutions of [START_REF] Szalai | Global dynamics of low immersion highspeed milling[END_REF] but not of [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF]. However, pathologically small choices will result in unnecessarily long computations.

Method of numerical investigation

To investigate the dynamics of the system, we perform oneparameter continuations in the chip width w, with all other parameters held constant, in order to compute the structure sketched in figure 1(b). In each computation, we begin with the trivial solution x(t) ≡ 0 and follow the branch of unstable periodic orbits that emanates from the subcritical Hopf bifurcation. These calculations are at first entirely smooth, since the oscillations represent small amplitude motions in which the tool does not leave the workpiece, and which can be followed [START_REF] Dombóvári | Estimates of the bistable region in metal cutting[END_REF] using the basic model [START_REF] Merdol | Multi frequency solution of chatter stability for low immersion milling[END_REF]. However, as the amplitude of the unstable periodic orbit grows, then it begins to incorporate motion in which the tool flies, for which the full model [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF] and its smoothing [START_REF] Szalai | Global dynamics of low immersion highspeed milling[END_REF] are required.

Note that linear stability analysis of the steady cutting solution x ≡ 0 leads to the well-known lobe structure in the Ω-w stability chart [1-3, 11, 15, 27], see figure 3. In this paper, we focus on computations which cut the 'first' lobe at its minimum with respect to w, for which the non-dimensional spindle speed Ω is highest, indicated by the vertical section in figure 3. (For the choice of physical parameters we have adopted in section 2, 5) which captures motions with no loss of contact, and the large circle on this branch denotes the switching point at which zero chip thickness is attained. In panel (a) we fix ε 2 = 10 -7 and reduce ε 1 , and in (b) we fix ε 1 = 5 • 10 -7 and reduce ε 2 . In these limits, the unstable portion of the branches accumulate on the branch for (5), and the fold points accumulate on the switching point for (5) if ε 1 is reduced sufficiently.

we thus have Ω 1.37.) In all plots we then present the chip width normalised by its minimum value w min at the Hopf bifurcation. Preliminary computations (not presented here) indicate similar bifurcation structure for other constant Ω sections. Some of this similarity can be explained by the standard DDE invariance property, which relates the lobes to each other. Specifically, if we have a periodic orbit φ(t) with period T for a particular time delay τ = τ 0 , then φ(t) is also a periodic orbit for time delay τ = τ 0 +nT where n = 1, 2, . . .. However, stability properties are not preserved by this transformation. In consequence, whereas the secondary bifurcations of the continuation through the first lobe include many period doubling bifurcations, the continuations through other lobes typically involve Neimark-Sacker bifurcations giving rise to quasi-periodic oscillations which cannot be continued further by DDE-BIFTOOL.

Effect of the smoothing parameters

We now follow the continuation procedure described in section 3.2, in order to investigate the dynamics of the smoothed system [START_REF] Szalai | Global dynamics of low immersion highspeed milling[END_REF]. Of particular interest is the choice of the smoothing parameters, ε 1 and ε 2 , and how small they must be chosen for the dynamics of the non-smooth model [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF] to be reproduced robustly. Furthermore, we wish to establish the nature of the transition to chattering dynamics, at point D in figure 1(b). These questions are the focus of figure 4, which summarises our results.

Figure 4 establishes two things. Firstly, as ε 1 and ε 2 are reduced towards zero, there is an accumulation in the branches computed from the smoothed model [START_REF] Szalai | Global dynamics of low immersion highspeed milling[END_REF]. In fact, this accumulation establishes only the self-consistency of the smoothed model, and not that it reproduces the dynamics of the full model (10) -a point that we will return to in section 4.

Secondly, figure 4 shows that the branch of unstable periodic orbits emanating at the Hopf bifurcation 'turns over' in a cyclic fold (Saddle Node of a Limit Cycle) at which it regains stability. The location of this fold may then be examined as ε 1 and ε 2 are successively reduced. When this limiting process is followed, we find that the fold for the smoothed model ( 12) tends to the point at which the orbits of the permanent contact model ( 5) attain loss of contact between the tool and the workpiece.

In consequence, it seems that we have established the nature of the point D in figure 1(b). The limiting behaviour of the smoothed model indicates that the branch of unstable periodic orbits of the full non-smooth model [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF] regains stability via a non-smooth cyclic fold, at which point its orbits simultaneously lose contact between the tool and the workpiece. This is a wellknown type of grazing bifurcation [START_REF] Bernardo | Piecewisesmooth dynamical systems: theory and applications[END_REF]Chapter 7]. However, as we shall now see, a detailed investigation of the secondary bifurcation structure beyond the fold point reveals a more complicated situation.

Details of chattering dynamics

We now extend our analysis to a more detailed study of chattering motions -that is the region denoted E in figure 1(b). There are two components in this work.

• Firstly, in section 4.1, we extend our numerical continuation of the smoothed model [START_REF] Szalai | Global dynamics of low immersion highspeed milling[END_REF] to consider secondary bifurcations that occur beyond the fold point that we have already examined. In summary, we discover a complicated bifurcation structure in which many branches of period doubled orbits also accumulate at the fold point, but 'from the right' in the sense of figure 4. This result indicates that the onset of chattering motion is not via a simple grazing cyclic fold, and that initial chattering motions are probably not simple periodic orbits.

• Thus secondly, in section 4.2, we perform time-domain simulations of the full non-smooth model [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF] in order to eliminate artefacts due to the smoothing process. In summary, we establish that there is an immediate 'jump' to chaotic dynamics, and we trace out how the dynamics evolve as the chip width is increased, demonstrating the correspondence between periodic windows in the time simulation and stable sections of periodic orbit in the numerical continuation.

Numerical continuation and secondary bifurcations

We now follow on from figure 4 to continue branches of periodic orbits to the smoothed model [START_REF] Szalai | Global dynamics of low immersion highspeed milling[END_REF], as the chip width w is increased beyond the non-smooth fold point. In particular, we focus on computing secondary bifurcations and the branches that emanate from them. Our experience with computations of the smoothed model indicates that qualitative dynamics are relatively robust to the choice of the smoothing parameter ε 2 , but the choice of ε 1 is critical. Hence throughout we fix ε 2 = 10 -7 and consider a sequence of computations as ε 1 is reduced. Figure 5 summarises the results.

Further continuation of the branches shown in figure 4 for ε 1 = 10 -7 demonstrates the opening up of a stable perioddoubled bubble, as seen in figure 5(a). However, as ε 1 is decreased further, see figure 5(b), a second stable period-doubled bubble emerges so as to introduce a window of unstable periodfour orbits. Moreover, the original left-hand period-doubling bifurcation has moved towards the fold. This process continues in figure 5(c) and (d) as ε 1 is reduced further and we begin to recover the converged bifurcation structure of the model. In short we discover an apparent cascade of period-doubling bifurcations, where the left-hand bifurcation of each new bubble approaches the fold point as ε 1 is successively reduced.

Furthermore, as ε 1 is decreased, the criticality of the perioddoubling bifurcations changes from supercritical to subcritical. The change in criticality introduces additional cyclic fold bifurcations into the branch structure, labelled A and B in figure 5. It appears that, in the limit of ε 1 → 0, these cyclic folds also become non-smooth fold points. However, no further bifurcating branches are found to emerge from these extra non-smooth fold points. The change in criticality is further highlighted in figure 6 which shows the change in the period-3 and period-6 branches as ε 1 is changed from 5 • 10 -9 to 10 -9 .

In consequence, it seems the onset of stable periodic chattering motion at the primary fold point, discovered in section 3.3, is an artefact of the smoothing process. Instead, it seems that the fold point is a more complicated type of grazing bifurcation from which infinitely many branches of unstable periodic orbits emanate. This is reminiscent of the 'Big Bang bifurcation' of Avrutin and Schanz [START_REF] V. Avrutin | On multi-parametric bifurcations in a scalar piecewise-linear map[END_REF] and as such we label this point B 3 . The emerging periodic orbits are then absorbed at period-doubling bifurcations at larger values of the chip width parameter. In fact, a more detailed investigation also reveals the emanation of branches of unstable period-three and period-six orbits which do not connect with the period-one/two/four structure (other than at the fold point) but which can be readily be found by numerical simulation and then continued.

This bifurcation structure, consisting of an infinite accumulation of unstable periodic orbits, has been discovered in sim-pler non-smooth dynamical systems [18, see examples in Chapter 7] and is associated with the instantaneous onset of localised chaotic dynamics. However, a recent parallel study [START_REF] Barton | Stability calculations for piecewise-smooth delay equations[END_REF] which works directly on the non-smooth model [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF], demonstrates that the period doubling cascade is an artefact of the smoothing, and that only the first two right-most period-doubling bifurcations, see figure 5(b), persist in the perfectly non-smooth model. The ε 1 , ε 2 → 0 limiting process thus has features which we could not explore fully from mathematical perspective, but the results obtained so far are also supported by the time-domain simulations of the limiting case.

Time-domain simulations

To establish the validity of the smoothed model ( 12) we resort to numerical simulation of the full non-smooth model [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF]. If the smoothed model is an accurate representation of the non-smooth model, then there should be good agreement between the results away from any non-smooth grazing bifurcations. Furthermore, we present the results for numerical simulations of the smoothed model [START_REF] Szalai | Global dynamics of low immersion highspeed milling[END_REF] to enable a complete comparison of the results.

Moreover, while numerical simulation cannot provide the full bifurcation structure of the model, it will reveal any chaotic attractors in the system which cannot be studied directly using numerical continuation methods. Thus numerical continuation and numerical simulation are complementary to each other in elucidating the dynamics of a model.

To integrate [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF] the Matlab routines dde23 and ode45 are used alternately to integrate the DDE and ODE parts of the system respectively; the built-in event detection features of the Matlab routines enable the switching points to be found to machine precision (since the switch is defined by an algebraic equation), thus the accuracy of the integration is unaffected by the switching. The algebraic constraint is implicitly incorporated into the system through modifications to the initial function segment provided to dde23. At each parameter value, [START_REF] Litak | Chaotic vibrations in a regenerative cutting process[END_REF] is integrated for 1000 non-dimensional time units and the first 200 time-units of data are then discarded to account for the transient behaviour. The output of each simulation run is used as the initial condition for the next simulation run.

Figure 8 shows results from numerical simulations of (10) overlaid with the cloud from numerical continuation of [START_REF] Szalai | Global dynamics of low immersion highspeed milling[END_REF]. As suggested by the numerical continuation results, the primary fold point / grazing bifurcation acts as an organising centre for the dynamics.

The values on the vertical axis of figure 8 are defined by

ξ := sgn(x(t)) x(t) 2 + x(t -τ) 2 , ( 13
)
at the point the orbit intersects the Poincaré section

Π := {x(t) + x(t -τ) = 0}. ( 14
)
This section is chosen for convenience as all periodic orbits pass through it (see figure 7). For the parameters chosen here, we have x(t -τ) = χ(t -τ) at the switching surface. Therefore, the switching surface Σ := {h(t) = 0} is orthogonal to Π in the -1 0 1

-1 0 1 x t ( ) orbit x t ( ) -τ S f S c S P x
Figure 7: Shows a two dimensional projection of phase-space and a representative period-one orbit. Π marks the Poincaré section defined by [START_REF] Kalmár-Nagy | Experimental and analytical investigation of the subcritical instability in metal cutting[END_REF] and Σ marks the switching surface that separates the cutting and flying vector-fields (S c and S f ).

(x(t), x(t-τ))-plane and it takes the constant value of ξ Σ = 1/ √ 2. Thus, grazing bifurcations (or approximations thereof) can be immediately seen as locations where the solution branches intersect (or in the smoothed case, come close to) ξ = ξ Σ .

At the primary grazing bifurcation (labelled B 3 in figure 8), from which many branches of periodic orbits emerge, a highly localised chaotic attractor also emerges. This chaotic attractor (cw a ) appears to be the global attractor of the system until the grazing bifurcation of the period-6 orbits (figure 8(b) where p 6 and p 3 were continued using ε = 10 -9 ). The numerical simulation also indicates that the period-3 branch continues further than the numerical continuation results show before terminating at a fold / grazing bifurcation (G 3 in figure 8(d)), thus extending the chaos-free window in parameter space. However, for the period-2 and period-1 branches the agreement between numerical continuation of the smoothed system and numerical simulation of the non-smooth system is excellent.

Figure 9 shows two dimensional projections of the Poincaré section defined by [START_REF] Kalmár-Nagy | Experimental and analytical investigation of the subcritical instability in metal cutting[END_REF]. Figures 9(a) and (c) show the phasespace when there are stable periodic orbits (solid symbols) and unstable periodic orbits (open symbols), and figures 9(b) and (d) show the phase-space when there is a stable chaotic attractor. In both cases, the agreement between simulation and continuation is good. In particular, all of the unstable periodic orbits can be seen to reside within the chaotic attractor as is expected.

Conclusion

This paper has been motivated by chattering behaviour in machining processes. We have focussed on a simplified one degree-of-freedom model of orthogonal cutting, as a representative example displaying the regenerative effect common to turning processes. Even this simplified modelling approach, when incorporating full chattering motions in which the tool loses contact with the workpiece, yields a non-smooth delay differential algebraic equation (DDAE) which cannot be analysed by standard numerical continuation tools.

Our chief result has been to derive a smoothed version of the orthogonal cutting system which is a pure delay differential equation (DDE) without algebraic effects. This is achieved by relaxing time to eliminate the algebraic effect, and by relaxing space to eliminate the non-smoothness that occurs at the loss of contact. These relaxations introduce two (small) smoothing parameters that govern the accuracy of the new approximate DDE system.

We have then analysed the smooth DDE system using the numerical continuation software DDE-BIFTOOL and we have compared the dynamics with initial value simulations of the full non-smooth DDAE. In broad qualitative terms, there is excellent agreement thus vindicating our approach. But the numerical continuation reveals dynamical structure and mechanisms which could not otherwise be investigated.

In particular, we have used numerical continuation to investigate the onset of chattering dynamics, at the point labelled D in figure 1(b). For the cubic cutting force characteristics used here, we have shown that the unstable branch of periodic orbits folds over at the same point in parameter space where its solutions first lose contact with the workpiece. This result has been suspected for some time (e.g., [START_REF] Dombóvári | Estimates of the bistable region in metal cutting[END_REF][START_REF] Stepan | Nonlinear Dynamics of Material Processing and Manufacturing[END_REF]) but the matter is now resolved. However, in computations not shown here, we have discovered alternative mechanisms for the onset of chattering when other cutting force relationships are used, such as the three-quarter rule. When the fine detail of solution branches is considered, some discrepancies can be found between the smoothed DDE and the full non-smooth DDAE system. Although we have established the convergence of the solution branches of the smoothed system in the limit of vanishing smoothing parameters, this converged structure is not always consistent with the DDAE in the neighbourhood of all non-smooth bifurcation points. Specifically, the route to chaos in the smoothed DDE is through a classical period-doubling cascade whereas recent parallel work has shown that the onset of chaos in the non-smooth DDAE is due to a grazing bifurcation [START_REF] Barton | Stability calculations for piecewise-smooth delay equations[END_REF]. Moreover, this grazing bifurcation (point D in figure 1(b)) acts as an organising centre from which all the chattering orbits emerge.

A further point to highlight is the apparent discrepancy between the results shown here and the results of Wahi and Chatterjee [START_REF] Wahi | Self-interrupted regenerative metal cutting in turning[END_REF] which arises due to the use of a different cutting force model. Here the cutting rule of Shi and Tobias [START_REF] Shi | Theory of finite-amplitude machine-tool instability[END_REF] is used whereas [START_REF] Wahi | Self-interrupted regenerative metal cutting in turning[END_REF] uses the three-quarter rule [START_REF] Kienzle | Spezifische schnittkräfte bei der metallbearbeitung[END_REF]. When using the three-quarter rule, no chaotic chattering motion is found, only periodic motion persists. As such, further experimental work is required to determine which of the force laws is more applicable. Indeed, the results presented here may provide a distinguishing feature which can be used to differentiate between the force laws, as direct force-displacement measurements cannot be performed with sufficient accuracy. We emphasize again, however, that experimental evidence in the literature [START_REF] Shi | Theory of finite-amplitude machine-tool instability[END_REF][START_REF] Stepan | Modelling nonlinear regenerative effects in metal cutting[END_REF] supports the existence of chaotic oscillations.

Finally, we should point out that our broad philosophy is fo-cussed on detailed parameter studies of simple low degree-offreedom models. The fine details that we discover may not survive real-world perturbations such as process damping at low spindle speeds [START_REF] Altintas | Identification of dynamic cutting force coefficients and chatter stability with process damping[END_REF] or the incorporation of higher order structural modes. However, given that very simple models display complex dynamics (including chaos) in region E, it seems plausible that the dynamics of more faithful models are at least as complicated -but their investigation remains for future work.

Figure 1 :

 1 Figure 1: Panel (a) shows a photograph of a machined workpiece. The smooth surface to the left was produced by steady cutting, whereas to the right we may observe so-called sunflower spirals which are a typical result of tool chatter (photograph provided by Ideko). Panel (b) depicts a schematic bifurcation diagram showing typical dependence of the tools motion on the chip width w. Here solid / dashed lines represent respectively stable / unstable steady-state cutting. From the Hopf point B an unstable period-one branch marked C emerges. Our interest is in point D, which denotes the onset of chattering motions at E.

Figure 2 :

 2 Figure 2: Panel (a) shows a schematic of the cutting tool in contact with the workpiece. The workpiece is rotating with angular velocity Ω and the tool enters the workpiece with feed velocity v f = h 0 Ω/(2π). Panel (b) is a graphical representation of the model (2) where the forces on the cutting tool are proportional to the chip thickness h(t) = q(t -τ)q(t) + h 0 , where τ = 2π/Ω is the period of revolution.

Figure 3 :

 3 Figure 3: Two-parameter bifurcation diagram of the steady-state cutting. The black curve represents a Hopf bifurcation curve; along this curve the bifurcations are subcritical. Consequently, there are unstable periodic orbits coexisting with the stable steady-state cutting solutions in the grey region of the diagram.

Figure 4 :

 4 Figure4: The panels show a series of one-parameter continuations (bold curves) of the smoothed model[START_REF] Szalai | Global dynamics of low immersion highspeed milling[END_REF] in the vicinity of the fold-type point at which the branches turn over. The grey curve indicates a continuation of the unstable periodic orbit of the basic model (5) which captures motions with no loss of contact, and the large circle on this branch denotes the switching point at which zero chip thickness is attained. In panel (a) we fix ε 2 = 10 -7 and reduce ε 1 , and in (b) we fix ε 1 = 5 • 10 -7 and reduce ε 2 . In these limits, the unstable portion of the branches accumulate on the branch for (5), and the fold points accumulate on the switching point for (5) if ε 1 is reduced sufficiently.

Figure 5 :

 5 Figure 5: Panels (a)-(d) show the convergence of the bifurcation structure of the smoothed model (12) as ε 1 → 0. Throughout we fix ε 2 = 10 -7 . In (a) there are two period-doubling bifurcations connected by a branch of period-two orbits. As ε 1 is decreased a 'bubble' of period-four orbits emerges as shown in (b). Panels (c) and (d) appear to show the converged bifurcation structure. All of the left-hand period-doubling bifurcations have moved towards the fold point, which suggests that all the branches shown actually emerge from the grazing bifurcation in the non-smooth model (cf. the labelled point D). Isolated branches of period-three and period-six orbits, connected by a period-doubling bifurcation at point C, are shown in (c) and (d) respectively.

Figure 6 :

 6 Figure 6: Panels (a) and (b) show a zoom of the marked region in figure 5( f ) for (a) ε 1 = 5 • 10 -9 and (b) ε 1 = 10 -9 . The period-doubling bifurcation connecting the period-3 and period-6 branches clearly changes criticality.

Figure 8 :

 8 Figure 8: Panel (a) and (c) show bifurcation diagrams calculated by numerical integration of the smoothed (12) and the non-smooth model (10), respectively. Panel (b) contains the branches of periodic orbits obtained by numerical continuation overlaid with the results of panel (c) for comparison. The smoothing parameter ε 1 was 5 • 10 -9 in panel (a) and (b) except for branches p 3 and p 6 in panel (b) where ε 1 was 10 -9. In all panels, the vertical axis is defined by (13) in combination with the Poincaré section[START_REF] Kalmár-Nagy | Experimental and analytical investigation of the subcritical instability in metal cutting[END_REF]. It can be seen that a localised chaotic attractor emerges from the primary fold point / grazing (Big Bang) bifurcation B3 . Panel (d) shows the relative locations of the different bifurcations, where PD denotes a period-doubling bifurcation and G denotes a grazing bifurcation.

Figure 9 :

 9 Figure 9: Two-dimensional projections of the Poincaré section (14) are shown in panels (a)-(d) corresponding to the parameter values marked in figure 8. These show the development of the bifurcation structure as the bifurcation parameter w is decreased. Panel (a) shows a stable period-two orbit (solid diamonds) and an unstable period-one orbit (open circle). Panel (b) shows a chaotic attractor with unstable period-one (open circle), period-two (open diamond) and period-four (open square) orbits. Panel (c) shows an example in the periodic window of figure 8 a stable period-three orbit (solid triangles), an unstable period-three orbit (open triangles) and the unstable period-one, -two and -four orbits. Panel (d) shows the chaotic attractor with the unstable periodic orbits embedded.

6. Acknowledgements Z.D. and G.S. were supported by the Hungarian Scientific Research Foundation OTKA Grant No. K68910, HAS-BUTE Research Group on Dynamics of Machines and Vehicles, Hungarian-Spanish Intergovernmental S&T Cooperation Programme, Grant No. OMFB-01265/2007 and the OPENAER project (CENIT program of CDTI). D.A.W.B was supported by a Great Western Research fellowship. R.E.W. was supported by EPSRC grant EP/E055567/1.