Tobias Torp-Pedersen

PhD Heather A Boyd

Gry Poulsen

Birgitte Haargaard

Jan Wohlfahrt

Jonathan M Holmes

Mads Melbye

Perinatal Risk Factors for Strabismus

Keywords: Strabismus, Reproductive factors, Birth weight, Gestational age, Congenital malformations, Head circumference

Background:

Little is known about the etiologic factors underlying strabismus. We undertook a large cohort study to investigate perinatal risk factors for strabismus, overall and by subtype.

Methods: Orthoptists reviewed ophthalmologic records for Danish National Birth Cohort children examined for strabismus in hospital ophthalmology departments or by ophthalmologists in private practice. Information on perinatal characteristics was obtained from national registers.

We used log-linear binomial regression and polytomous logistic regression to estimate risk ratios for strabismus overall and by strabismus subtype, respectively.

Results: Among 96 842 DNBC children born in Denmark between 1996 and 2003, we identified 8783 children who had been evaluated for strabismus. Ophthalmologic records were available for 5655 of these children, of whom 1321 were diagnosed with strabismus. In multivariable analysis, low birth weight, prematurity, large head circumference and presence of congenital abnormalities were all associated with increased risk of strabismus. Presence of congenital abnormalities was more strongly associated with exotropia than with esotropia. Of 183 exotropia cases, 40 (22%) had a congenital abnormality. While not associated with esotropia, delivery by Cesarean section was associated with exotropia (RR=1.65; 95% CI 1.16-2.34). After adjustment for birth weight, Apgar score at five minutes, multiple gestation and parental ages were not associated with strabismus overall.

Conclusions: Congenital abnormalities, low birth weight, prematurity, and large head circumference were independent risk factors for strabismus. Differences in risk factors for esotropia and exotropia suggest that strabismus subtypes may have different underlying etiologies. The proportion of exotropic children with congenital abnormalities suggests that a large angle constant exotropia in an infant should alert physicians to the possibility of a congenital malformation.

Strabismus is a common disorder of the eyes, with a reported childhood prevalence of 2%-6% in Western countries. [START_REF] Mcneil | Patterns on visual defects in children[END_REF][START_REF] Frandsen | Some results from a clinical-statistical survey on strabismus among Copenhagen children[END_REF][START_REF] Chew | Risk factors for esotropia and exotropia[END_REF][START_REF] Nordlöw | Squint -The frequency of onset at differenet ages, and the incidence of some associated defects in a Swedish population[END_REF][START_REF] Graham | Epidemiology of strabismus[END_REF][START_REF] Williams | Prevalence and risk factors for common vision problems in children: data from the ALSPAC study[END_REF] Left uncorrected, strabismus can result in amblyopia, which may lead to permanent loss of vision.

Little is known about the etiologic factors underlying strabismus, which probably encompasses multiple etiologically distinct entities. Many types of strabismus are typically detected in childhood, suggesting that events occurring during the perinatal period might be important.

Several investigators have reported low birth weight as a risk factor for strabismus. (3;7-14) Although some studies have found prematurity to be associated with strabismus, (3;6;11;13-15) no study has demonstrated an association between gestational age at birth and strabismus after adjusting for birth weight. (14;16;17) Studies addressing other perinatal risk factors for strabismus are limited. (3;6;13;17) The sparse scientific literature on perinatal risk factors for strabismus prompted us to undertake a large population-based cohort study to address this issue. Specifically, we investigated the effect of maternal and paternal age, maternal pre-eclampsia, multiple gestations, birth weight, gestational age, head circumference, Cesarean delivery, neonatal asphyxia, Apgar score 5 minutes post-delivery, neonatal septicemia, and the presence of congenital abnormalities on the risk of strabismus, overall and by strabismus subtype.

Methods

Data Sources

Our study cohort consisted of the 96 842 children born alive between 1996 and 2003 to women participating in the Danish National Birth Cohort (DNBC). (18;19) We collected information from national registers concerning variables of interest and the whereabouts of ophthalmologic records, if any, for each child. We used the unique personal identification number (PIN),

Medical Birth Register

The Medical Birth Register has recorded live births and stillbirths in Denmark since 1973. The register contains information on birth characteristics (e.g. birth weight, gestational age, Apgar score at 5 minutes), birth outcomes, and delivery complications, recorded by midwives attending the delivery.

National Patient Register

Reporting of data on Danish residents' hospital admissions and outpatient visits has been mandatory since 1977 and 1994, respectively. The National Patient Register contains detailed information on date of hospitalization or outpatient visit, codes for discharge diagnoses, and codes for any surgical procedures. [START_REF] Andersen | The Danish National Hospital Register. A valuable source of data for modern health sciences[END_REF] Since 1994, the National Patient Register has coded discharge diagnoses using the 10 th Revision of the International Classification of Diseases (ICD-10) and surgical procedures using Nordic Medico-Statistical Committee Classification of Surgical Procedures (NCSP) codes.

Health Security System

Health care is free of charge for all residents of Denmark. The Danish health care system reimburses specialists in private practice for their services. The Health Security System records these reimbursements and contains detailed information on procedures and examinations performed (although no information on specific diagnoses). onset before the age of 6 months that could not be corrected with glasses. Accommodative esotropia was defined as an esotropia that could be corrected with glasses, either completely (fully accommodative) or partially (partially accommodative). In our analyses, fully and partially accommodative esotropia were combined as "accommodative esotropia". We classified esotropia that did not meet the criteria for congenital or accommodative esotropia as "other esotropia". "Exotropia" included all exotropia diagnoses. We classified children with any orbital nerve paresis as having "other strabismus".

Case Ascertainment

Ascertainment of Perinatal Characteristics

Using the Medical Birth Register, we obtained information on mode of delivery, birth weight, gestational age at birth, head circumference, Apgar score at 5 minutes, single/multiple gestation, neonatal septicemia, neonatal asphyxia, maternal pre-eclampsia, and parental ages. We searched the National Patient Register for ICD-10 codes for asphyxia (I11-69, J, O68, P20-29, R09) and septicemia (P36) in the child up to 20 days after delivery. We also searched for congenital abnormality diagnoses (ICD-10 codes Q00-Q99) at any age. We grouped congenital abnormalities based on proximity to the eyes, ranked as follows: 1) chromosomal abnormalities;

2) malformations of the head and neck (not including cleft palate/lip); 3) other malformations (including cleft palate/lip); and 4) no malformations. Children with multiple abnormalities were classified based on the diagnosis ranking highest in the above list.

Data Analysis

We used log-linear binomial regression to estimate risk ratios (referred to hereafter as relative risks) associated with individual perinatal risk factors for strabismus overall. All estimates were adjusted for birth year to take into account differences in follow-up time according to birth year.

We also examined the effect of adjusting for birth weight and restricting our analyses to children without congenital abnormalities. In addition to this, we then performed a multivariable analysis including only children without congenital abnormalities; to be conservative, we initially included in this analysis birth weight, birth year, and variables with P<0.10 in the individual variable analyses after adjustment for birth year and birth weight and exclusion of children with congenital abnormalities. Using a backward elimination procedure, we then reduced the multivariable model to include only those variables that were independent risk factors (P<0.05) for strabismus overall.

We performed additional multivariable analyses adjusting for the highest social class of the household (academics/employed requiring highly specialized skills, employed requiring medium skills, skilled labor, students, unskilled labor, and unemployed), daily maternal smoking dose (0, 0.1-4.9, 5-9.9, 10+ cigarettes per day) and weekly maternal alcohol consumption (0, 0.1-0.9, 1-2.9, 3-4.9, 5+ 12-gram units of alcohol per week) during pregnancy. This was done primarily to evaluate whether these risk factors were confounders of the observed perinatal risk factorstrabismus associations,(21) but also to evaluate whether confounding was introduced by adjusting for birth weight in analyses where birth weight could be regarded as an intermediate variable. [START_REF] Cole | Fallibility in estimating direct effects[END_REF] We used polytomous logistic regression to compare risk factors for strabismus subtypes. To ensure that we did not overlook variables of importance to specific subtypes, risk factors with P<0.2 in the analyses of strabismus overall with adjustment for birth year and birth weight and exclusion of congenital abnormalities were considered in the comparison of strabismus subtypes. We decided a priori to adjust strabismus subtype analyses for birth weight.

We used SAS Version 9.2 (SAS Institute, Inc., Cary, NC, USA) for data analysis.

Results

F

o r R e v i e w O n l y 9
Evaluation of 5655 ophthalmologic records yielded 1321 strabismus cases. One child with strabismus caused by trauma was excluded, leaving 1320 cases for analysis in a cohort of 96 842 children. We first present associations between perinatal risk factors and risk of strabismus overall, followed by the results of separate analyses for congenital esotropia, accommodative esotropia, all exotropia, and other strabismus.

Figure 1a shows the association between birth weight and any strabismus, adjusted only for birth year. Children with birth weights below 2500g had an increased risk of strabismus, compared with children with birth weights above 2500g (P<0.001). Additional adjustment for gestational age affected the risk estimates for birth weight very little (P=0.001) (Figure 1b).

Figure 1c shows the association between gestational age and any strabismus, adjusted only for birth year; children born before week 37 of gestation had an appreciably higher risk of strabismus than children born in weeks 37-41 (P<0.001). However, after adjusting for birth weight, the association with gestational age diminished markedly (P=0.05) (Figure 1d) Table 1 presents associations between perinatal characteristics and strabismus overall. Relative risk estimates are adjusted for 1) birth year only, 2) birth year and birth weight, and 3) birth year and birth weight in children with no malformations or chromosomal abnormalities. When only adjusting for birth year, asphyxia, a low Apgar score, Cesarean delivery, small and large head circumferences, maternal pre-eclampsia, multiple gestation, and neonatal septicemia were all associated with strabismus. We obtained similar results when we further adjusted the estimates for social class and maternal smoking and alcohol consumption during pregnancy (data not shown). After additional adjustment for birth weight, children with neonatal asphyxia remained at elevated risk of strabismus (Table 1). Similarly, children born by Cesarean section had an increased risk of strabismus compared with children born vaginally (Table 1). Estimate magnitudes for neonatal asphyxia and Cesarean section remained similar when children with congenital abnormalities were excluded from the analyses (Table 1).

Compared with children born with a head circumference of 34-37 cm, children with larger head circumferences had a 36% (95% CI 11%-69%) higher risk of strabismus, after additional adjustment for birth weight. The risks associated with a head circumference of 33 cm or less remained elevated, although more modestly so, after adjustment for birth weight. Our head circumference estimates did not change further when children with congenital abnormalities were excluded from the analyses (Table 1).

To explore the robustness of the above results, we performed an additional multivariable analysis (with backward elimination of variables), initially including the variables Cesarean delivery, maternal pre-eclampsia, head circumference, gestational age and birth weight, in children without congenital abnormalities (results not shown). Using this approach, we also found that only birth weight (P homogeneity =0.001), gestational age (P homogeneity =0.05), and head circumference (P homogeneity =0.02) were associated with strabismus risk. Compared with children weighing 3000-3500 g at birth, children weighing less than 2000g had a 42% (95% CI -15%-135%) increased risk of strabismus, while children weighing 2000-2500g had a 74% (95% CI 28%-137%) increased risk of strabismus. Compared with children born in gestational weeks 37-42, children born in week 32 or earlier had a 19% (95% CI -34%-114%) increased risk of strabismus, while children born in weeks 33-36 had a 39% (95% CI 8%-80%) increased risk of strabismus. Children born post-term had an 18% (95% CI -2%-34%) decreased risk of strabismus. 2 presents associations between congenital abnormalities and risk of strabismus. After adjustment for birth weight, children with abnormalities had a 57%-655% increased risk of strabismus compared with children with no abnormalities, depending on the abnormality. Table 3 presents associations between perinatal risk factors and esotropia (n=1048) and exotropia (n=183) separately. All estimates are adjusted for birth weight. Congenital abnormalities of all types were more strongly associated with exotropia than esotropia (P<0.001). In fact, 40 (22%) of the 183 exotropia cases also had one or more congenital abnormalities, compared with 12% of esotropia cases. The most frequent abnormalities among the 30 exotropic children with malformations in only a single organ system were musculoskeletal (n=9), facial (n=7), and malformations of the nervous system (n=4).

Compared with children delivered vaginally, children born by Cesarean section had no increased risk of esotropia, while their risk of exotropia was increased by 65% (95% CI 16%-134%).

Maternal pre-eclampsia, neonatal septicemia, gestational age at birth, and head circumference did not affect the risks of esotropia and exotropia differently (Table 3). We could not obtain ophthalmologic records for 3128 children. However, these children did not differ from children with available records with respect to the distribution of birth year, birth weight, or gestational age at birth (data not shown), and participating and non-participating ophthalmologists did not differ in terms of the numbers of records we were attempting to access or the region of the country in which they practiced (data not shown). Children with missing records were excluded from our primary analyses. However, we also performed an alternate analysis where these children were considered strabismus-free. This had very little effect on the effect estimates and did not change our conclusions (data not shown).

Discussion

In this population-based cohort study, strabismus risk was associated with low birth weight, low gestational age at birth, large head circumference and congenital abnormalities. Congenital abnormalities were more strongly associated with exotropia than with esotropia; of particular note, as many as 22% of exotropic children also had a malformation or chromosomal abnormality. Children born by Cesarean section had an increased risk of exotropia, while their risk of esotropia was unaffected.

The inverse association we observed between strabismus and birth weight is in line with previous findings in a cohort of children born in 1959-1966.(3) Other studies have reported strabismus prevalences of 12%-36% in very low birth weight children (birth weight less than 1500 g),(7;9;10;12) compared with 2%-6% in the general pediatric population.(1-6) The Sydney Myopia Study, a large population-based cross-sectional study of ophthalmologic conditions in young children (n=1739, 48 with strabismus when assessed at 6-7 years of age), found that children weighing less than 2500 g at birth were 3.5 times more likely than children who were greater than 2500 g at birth to have any type of strabismus; [START_REF] Robaei | Factors associated with childhood strabismus[END_REF] 14) Compared with children of normal birth weight, we found a 2-fold increased risk of strabismus both in children weighing less than 2000 g and in those weighing 2000-2500 g at birth. The large size of our study cohort resulted in robust estimates, with an upper confidence limit of less than 2.6 for both estimates. In contrast to our findings that low birth weight (< 2500 g) increased the risk of esotropia slightly more than 2-fold and the risk of exotropia only modestly (approximately 20-50%), the Sydney Myopia Study found that low birth weight children were 3.4 times more likely to be esotropic and 2.5 times more likely to be exotropic than children with normal birth weights, although both estimates were based on very small numbers of children.( 13) Major and colleagues reported an unadjusted odds ratio for congenital esotropia of 5.9 for children weighing less than 2500 g,( 11) which compares favorably with our relative risks of 2.66 and 3.58 for children weighing less than 2000 g and 2000-2499 g, respectively.

Several groups have reported an association between prematurity and an increased risk of strabismus,(3;11;13-15) but only one group adjusted their gestational age estimate for birth weight, reporting a residual effect of gestational age after adjustment.(3) In their report on vision problems in 7538 7-year-old children from the Avon Longitudinal Study of Parents and Children (ALSPAC), Williams and colleagues examined the effects of gestational age alone and birth weight adjusted for gestational age.(6) They found that prematurity (gestation < 37 weeks) increased the risk of "clinically significant convergent strabismus" (esotropia) by roughly 2.5fold, but did not affect the risk of "clinically significant divergent strabismus" (exotropia), whereas intrauterine growth retardation (birth weight adjusted for gestational age more than 2 standard deviations lower than mean) increased the risk of exotropia 4.3-4.5 times but had no effect on the risk of esotropia. [START_REF] Williams | Prevalence and risk factors for common vision problems in children: data from the ALSPAC study[END_REF] In our analyses, the influence of birth weight appeared to explain most of the effect of gestational age, whereas adjusting for gestational age did not change the association between birth weight and strabismus. We found very strong associations between several groups of congenital abnormalities and strabismus, whereas previous studies have only found associations between specific abnormalities and strabismus. [START_REF] Creavin | Ophthalmic abnormalities in children with Down syndrome[END_REF][START_REF] Casteels | Ocular findings in children with a microdeletion in chromosome 22q11.2[END_REF][START_REF] Kuchukhidze | Hypoplasia of deep cerebellar nuclei in joubert syndrome[END_REF][START_REF] Mansour | Ocular pathology in congenital heart disease[END_REF][START_REF] Izquierdo | Strabismus in the Marfan syndrome[END_REF] In our study, congenital abnormalities were more strongly associated with exotropia than with esotropia, with 22% of exotropic children also having one or more abnormalities. Similarly, Izquierdo and colleagues reported a 5.6:1 exotropia:esotropia ratio in 573 subjects with Marfan syndrome,( 27) although esotropia is more common than exotropia in Western countries. (3;6;28) and Hunter and Ellis reported a higher association of genetic disorders particularly with infantile constant exotropia. [START_REF] Hunter | Prevalence of Systemic and Ocular Disease in Infantile Exotropia[END_REF] The large difference in the effect of congenital abnormalities on exotropia versus esotropia risk found in our study suggests that the pathogenesis of exotropia differs from that of esotropia. In support of this view, we also found that Cesarean delivery affected the risk of strabismus subtypes differently, increasing the risk of exotropia but not esotropia.

We found an association between large head circumference and strabismus risk, even after excluding children with congenital malformations. A large head circumference can point to both chromosomal abnormalities and conditions such as hydrocephalus, which could explain our findings. We found a similar association with small head circumference. Small neonatal head circumference has previously been described as a predictor of poor mental and motor development [START_REF] Badr | Predictors of neurodevelopmental outcome for preterm infants with brain injury: MRI, medical and environmental factors[END_REF] and may at least in part explain the observed association. Danish general practitioners were responsible for recruiting pregnant women for the DNBC; 64% of general practitioners nationwide participated in recruitment activities. More than 60% of eligible pregnant women receiving prenatal care from participating general practitioners enrolled in the DNBC, with the result that the DNBC included roughly 30% of all children born in Denmark during the enrollment period. Since neither a physician's decision to recruit women for the DNBC nor a pregnant woman's decision to participate in the cohort were likely to have been related to subsequent birth characteristics, our results are unlikely to have been subject to selection bias. Furthermore, a 2006 study concluded that although DNBC women were marginally healthier than the background population, any bias of risk estimates due to a "healthy woman" phenomenon would be small. [START_REF] Nohr | Does low participation in cohort studies induce bias?[END_REF] We were unable to retrieve ophthalmologic records for 36% of children who had been seen by an ophthalmologist and could have been diagnosed with strabismus. However, children with missing records were similar to those whose records we evaluated with respect to birth year, birth weight, and gestational age at birth. Furthermore, missingness was almost exclusively explained by OPP refusal to participate or by our inability to reach certain OPPs. Given the structure of the Danish health care system, we find it difficult to identify a mechanism whereby such non-participation could seriously bias our results.

Our study had several strengths. It was based on a large population-based cohort, which gave us the ability to study strabismus subtypes. Review of ophthalmologic records was performed by personnel specifically trained in strabismus classification, providing an accurate classification of strabismus. The use of national registers to which reporting is mandatory helped to ensure complete strabismus case ascertainment. Furthermore, health care in Denmark is free, allowing all residents equal access to medical care and ensuring that most strabismus cases get diagnosed.

Information on potential perinatal risk factors was registered in national registers by attending midwifes prior to, and independent, of strabismus diagnosis, minimizing potential recall bias.

To conclude, in a large population-based cohort study, we identified low birth weight, prematurity, large head circumference, and congenital abnormalities as risk factors for strabismus overall. Both congenital abnormalities and Cesarean delivery influenced the risk of Birth weight (g) RR q q q q q q q 0 1000 2000 3000 4000 5000

3 Medical

 3 Subject Headings: strabismus, risk factors, birth weight, prematurity, Cesarean section, congenital abnormalities, head circumference

  Danish residents (Danish citizens and persons with Danish residence permits) by the Danish Civil Registration System since 1968, to link individual-level information from the national registers for each cohort child.

  information with information from the National Patient Register (information available up to March 16 th , 2006) and the Health Security System (information available up to December 31 st , 2005) to identify DNBC children diagnosed with or suspected of having strabismus. Specifically, we searched the National Patient Register for strabismus diagnoses (ICD-10 codes H49-H51) and strabismus operations (NCSP codes KCEA-KCEW) to identify children with a hospital-confirmed strabismus diagnosis and children who had undergone strabismus surgery. We also searched the Health Security System for the strabismus evaluation code (19.2001) to identify children who had been evaluated for strabismus by an ophthalmologist in private practice (OPP).The search of the National Patient Register yielded 447 children with strabismus-related visits to ophthalmology and/or pediatric departments. The search of the Health Security System identified 8674 children who had visited 186 different OPPs. We subsequently requested ophthalmologic records for all children identified by the two searches. The number of children for whom records were requested (n=8783) does not equal the sum of the children identified in the two searches (9121) because there was some overlap between children identified using the two methods. We obtained records for 5655 (64%) of 8783 children evaluated for strabismus by an OPP and/or diagnosed with strabismus at a hospital department. We could not obtain records for 3127 seen by 66 OPPs, due mainly to OPP refusal to participate (n=31), non-response (n=23), or retirement (n=10). Only one of the missing records was a hospital chart.Two experienced orthoptists evaluated the ophthalmologic records and recorded specific strabismus subtype diagnosis (if any), date of diagnosis, and any concomitant disease relevant to strabismus (central nervous system, neurologic or facial bone disease). The orthoptists classified strabismus into the following groups: exotropia, congenital esotropia, accommodative esotropia, other esotropia, and other strabismus. Congenital esotropia was defined as an esotropia with

Table 4

 4 presents associations between perinatal risk factors and esotropia subtypes. Compared with children with no congenital abnormalities, children with chromosomal abnormalities or syndromes had an elevated risk of all three subtypes of esotropia. Children with malformations of the head and neck and children with other malformations had elevated risks of congenital esotropia and other esotropia, but not of accommodative esotropia. Birth weight was the only other variable that affected esotropia risk differently depending on subtype (p=0.04), with congenital esotropia yielding the highest relative risk estimates for low birth weight.

  modest low birth weight (1500-2499 g) increased the risk of any strabismus at 6-7 years 3.7-fold, compared with normal birth

  differently, with congenital abnormalities greatly increasing the risk of exotropia while only modestly increasing the risk of esotropia, and Cesarean delivery only affecting the risk of exotropia. These differences suggest that exotropia and esotropia have different underlying etiologies. A large angle constant exotropia in an infant should alert the physician to the possibility of an undetected congenital malformation.

Figure 1 .

 1 Figure 1.

Table 1 .

 1 Relative risk of strabismus according to birth characteristics in a cohort of 96,842 Danish children born between 1996 and 2003.

	Relative risk (95% confidence interval)
	Adjusted for birth weight and birth year

Table 2 .

 2 Relative risk of strabismus according to type of congenital abnormality in a cohort of 96 842 Danish children born between 1996 and 2003

	Relative risk (95% confidence interval)

Table 3 .

 3 Relative risks of esotropia and exotropia according to birth characteristics in a cohort of 96 842 Danish children born between 1996 and 2003

			All esotropia	All exotropia	Test for
							difference
							between
	Birth characteristics		Number of cases (n=1048)	Relative risk* (95% confidence interval)	Number of cases (n=183)	Relative risk* (95% confidence interval)	relative risks for esotropia and exotropia
	Birth weight	< 2000	41	2.20 (1.60-3.05)	5	1.52 (0.61-3.81)
	(g)	2000-2499	65	2.35 (1.80-3.07)	6	1.23 (0.53-2.86)
		2500-2999	126	1.29 (1.04-1.58)	22	1.26 (0.77-2.08)
		3000-3499	287	1 (ref)	51	1 (ref)
		3500-3999	296	0.85 (0.72-0.99)	59	0.95 (0.65-1.38)	P=0.73
		4000-4499	171	0.97 (0.80-1.17)	30	0.95 (0.61-1.49)
		4500+	54	1.21 (0.91-1.61)	9	1.12 (0.55-2.28)
		Missing	8	-	1	-
		P homogeneity		P<0.001		P=0.875
	Cesarean	No	840	1 (ref)	138	1 (ref)
	delivery	Yes	208	1.06 (0.91-1.25)	45	1.65 (1.16-2.34)	P=0.03
		P homogeneity		P=0.44		P=0.01
	Congenital	None	928	1 (ref)	143	1 (ref)
	abnormalities	Other	81	1.34 (1.06-1.68)	20	2.24 (1.40-3.57)
	Head / neck	21	2.80 (1.83-4.28)	13	11.7 (6.67-20.7)	P<0.001
	Chrom / syn **	18	5.47 (3.48-8.60)	7	16.1 (7.51-34.6)
		P homogeneity		P<0.001		P<0.001
	Gestational age	< 33	28	1.21 (0.69-2.12)	1	0.13 (0.01-1.17)
	at birth	33-36.9	100	1.39 (1.07-1.81)	10	0.70 (0.32-1.53)
	(weeks)	37-41.9	841	1 (ref)	156	1 (ref)	P=0.16
		≥ 42	79	0.81 (0.65-1.03)	16	0.90 (0.53-1.52)
		P homogeneity		P=0.03		P=0.19
	Head	≤ 30	27	1.41 (0.81-2.46)	1	0.13 (0.01-1.27)
	circumference	31-33	172	1.16 (0.95-1.41)	23	0.84 (0.50-1.42)
	(cm)	34-37 ≥ 38	696 100	1 (ref) 1.42 (1.13-1.79)	137 15	1 (ref) 1.08 (0.60-1.93)	P=0.17
		Missing	53	-	7	-
		P homogeneity		P=0.01		P=0.24
	Maternal pre-	No	991	1 (ref)	171	1 (ref)
	eclampsia	Yes	57	1.17 (0.89-1.54)	12	1.73 (0.95-3.14)	P=0.25
		P homogeneity		P=0.28		P=0.10
	Neonatal	No	1019	1 (ref)	179	1 (ref)
	septicemia	Yes	29	1.48 (1.02-2.15)	4	1.29 (0.47-3.54)	P=0.81
		P homogeneity		P=0.05		P=0.63
	* Adjusted for birth year and birth weight
	** Chromosomal aberrations and syndromes

Table 4 .

 4 Relative risks for esotropia subtypes according to birth characteristics in a cohort of 96 842 Danish children born between 1996 and 2003.

			Congenital esotropia	Accommodative esotropia	Other esotropia	Test for
									difference
	Birth characteristics	Number of cases (n=217)	Relative risk* (95% confidence interval)	Number of cases (n=432)	Relative risk* (95% confidence interval)	Number of cases (n=399)	Relative risk* (95% confidence interval)	between esotropia subtype relative risks
		< 2000	10	2.66 (1.36-5.20)	16	2.06 (1.22-3.45)	15	2.12 (1.24-3.63)
		2000-2499	20	3.58 (2.16-5.95)	17	1.47 (0.89-2.45)	28	2.66 (1.76-4.02)
		2500-2999	30	1.51 (0.97-2.35)	57	1.39 (1.02-1.91)	39	1.05 (0.73-1.51)
	Birth weight (g)	3000-3499 4000-4499 3500-3999	47 11 40	1 (ref) 1.12 (0.75-1.68) 0.67 (0.45-1.68)	115 23 80	1 (ref) 1.08 (0.82-1.44) 0.79 (0.61-1.01)	134 20 51	1 (ref) 0.76 (0.55-1.06) 1.01 (0.78-1.30)	P=0.04
		4500 +	58	1.22 (0.64-2.32)	120	1.23 (0.79-1.92)	109	1.18 (0.73-1.90)
		Missing	1	-	4	-	3	-
		p homogeneity		P<0.001		P=0.001		P<0.001
	Cesarean delivery	No p homogeneity Yes	170 47	1 (ref) P=0.65 1.08 (0.77-1.53)	348 84	1 (ref) P=0.46 1.10 (0.86-1.41)	322 77	1 (ref) P=0.90 1.02 (0.78-1.33)	P=0.91
		None	186	1 (ref)	400	1 (ref)	342	1 (ref)
		Other	20	1.56 (0.98-2.48)	19	0.72 (0.46-1.15)	42	1.81 (1.31-2.51)
	Congenital abnormalities	Head/neck Chrom/syn *	5 6	3.05 (1.26-7.41) 7.43 (3.29-16.8)	3 10	0.91 (0.29-2.82) 6.80 (3.63-12.7)	13 2	4.57 (2.64-7.93) 1.56 (0.39-6.24)	P=0.002
		*		P<0.001		P<0.001		P<0.001
		p homogeneity
		< 33	8	1.44 (0.48-4.29)	7	0.64 (0.24-1.72)	13	1.88 (0.77-4.55)
	Gestational age at birth (weeks)	33-36.9 ≥ 42 37-41.9	23 14 172	1.11 (0.64-1.95) 0.71 (0.41-1.23) 1 (ref)	39 34 352	1.44 (0.95-2.19) 0.82 (0.58-1.18) 1 (ref)	38 31 317	1.54 (1.00-2.37) 0.86 (0.59-1.25) 1 (ref)	P=0.64
		p homogeneity		P=0.55		P=0.08		P=0.18
		≤ 30	6	1.41 (0.45-4.42)	12	1.45 (0.61-3.47)	9	1.37 (0.54-3.49)
	Head circumference (cm)	31-33 Missing 34-37 ≥ 38	43 14 132 22	1.28 (0.84-1.96) -1 (ref) 1.58 (0.95-2.64)	71 15 293 41	1.25 (0.93-1.70) -1 (ref) 1.30 (0.91-1.87)	58 24 271 37	0.99 (0.70-1.38) -1 (ref) 1.47 (1.00-2.16)	P=0.93
		p homogeneity		P=0.02		P=0.25		P=0.25
	Maternal pre-eclampsia	No p homogeneity Yes	205 12	1 (ref) P=0.70 1.13 (0.62-2.04)	406 26	1 (ref) P=0.19 1.34 (0.88-2.04)	380 19	1 (ref) P=0.95 1.01 (0.63-1.64)	P= 0.24
	Neonatal septicemia	No p homogeneity Yes	208 9	1 (ref) P=0.05 2.12 (1.07-4.20)	422 10	1 (ref) P=0.49 1.26 (0.67-2.39)	389 10	1 (ref) P=0.39 1.34 (0.71-2.54)	P=0.51

* Adjusted for birth year and birth weight ** Chromosomal aberrations and syndromes

Acknowledgements

Financial support: The study was supported by grants from the Danish Eye Association [12-14-2006 ]; the Danish Medical Research Council [09-060139 ]; the Velux Foundation [05- ]; the Margrethe and Johs. F la Cour Foundation [04- ; and the Dagmar Marshall Foundation [11-03-2006 ].

Key Messages

In this population-based cohort study, low birth weight, prematurity, large head circumference and congenital abnormalities were associated with risk of strabismus overall.

Congenital abnormalities were much more strongly associated with exotropia than with esotropia. Given the substantial proportion of exotropic children with congenital abnormalities, a large angle exotropia in an infant should alert physicians to the possibility of a genetic disorder.

Children born by Cesarean section had an increased risk of exotropia, while their risk of esotropia was unaffected.

The differences in the effects of congenital abnormalities and Cesarean delivery on esotropia and exotropia suggest that strabismus subtypes have different underlying etiologies.


Andrew Yu Morozov

Emergence of Holling type III zooplankton functional response: Bringing together field evidence and mathematical modelling

Food-web population models are rather sensitive to parameterization of functional response in predation terms. Theoretical studies would predict enhancing of ecosystems' stability for a functional response of sigmoid type (Holling type III). The choice of a correct type of response is especially important for modelling outcome of grazing control of algal blooms by zooplankton in nutrient-rich ecosystems. Extensive experiments on zooplankton feeding in laboratories show non-sigmoid nature of response for most herbivorous zooplankton species.

As a consequence, there is a strong opinion in the literature that the implementation of Holling III type grazing in plankton models is biologically meaningless. I argue, however, that such an 'evident' claim might be wrong and sigmoid functional responses in real plankton communities would emerge more often than it was suggested earlier. Especially, this concerns plankton models without vertical resolution, which ignore heterogeneity in species vertical distribution. By having conducted extensive literature search of data on zooplankton feeding in situ, I show that vertical heterogeneity in food distribution as well as active food searching behaviour of zooplankton can modify the type of functional response.

Introduction

Dynamics of food-web population models largely depends on parameterization of functional response in predation terms (Oaten and Murdoch, 1975a,b;[START_REF] Hassell | Sigmoid functional responses and population stability[END_REF][START_REF] Murdoch | Plankton abundance and dynamics across nutrient levels: tests of hypotheses[END_REF][START_REF] Gross | Enrichment and food-chain stability: the impact of different forms of predator-prey interaction[END_REF]. This is especially true for modelling of eutrophic aquatic ecosystems where the choice of type of zooplankton functional response becomes crucial for predicting the outcome of grazing control of algae [START_REF] Truscott | Ocean Plankton Populations as Excitable Media[END_REF][START_REF] Hernández-Garcõa | Sustained plankton blooms under open chaotic flows[END_REF][START_REF] Fussmann | Community response to enrichment is highly sensitive to model structure[END_REF]. Implementation of an appropriate type of functional response for herbivorous zooplankton and, in particular, the correctness of the use of Holling type III (a sigmoid-shaped) response has provoked a vivid discussion in the literature [START_REF] Rothhaupt | Growth-rate de-pendent feeding rates in Daphnia pulicaria and Brachionus rubens: adaptation to intermediate time-scale variation in food abundance[END_REF]Scheffer and De Boer, 1996;[START_REF] Murdoch | Plankton abundance and dynamics across nutrient levels: tests of hypotheses[END_REF][START_REF] Jeschke | Consumer-food systems: why type I functional responses are exclusive to filter feeders[END_REF][START_REF] Sarnelle | Type III functional response in Daphnia[END_REF]. In this paper, I will consider the current problem from a different angle.

The zooplankton functional response is usually defined as the specific food intake rate (i.e.

per zooplankton biomass, per unit of time) as a function of the ambient food density.

Conventionally, there exist three main types of functional response known as Holling types

A c c e p t e d m a n u s c r i p t 3 [START_REF] Jeschke | Predator functional responses: discriminating between handling and digesting prey[END_REF][START_REF] Gentleman | Functional responses for zooplankton feeding on multiple resources: a review of assumptions and biological dynamics[END_REF]. These are Holling type I (linear response), Holling type II (concave upward response) and Holling type III (sigmoid response having a concave downward part at low food density). Note that replacing a functional response in a model with another one of a different type would completely modify models' properties (e.g. [START_REF] Gross | Enrichment and food-chain stability: the impact of different forms of predator-prey interaction[END_REF][START_REF] Fussmann | Community response to enrichment is highly sensitive to model structure[END_REF]. For instance, some models of plankton blooms simply do not 'work' properly in case of replacing Holling type III with a non-sigmoid one (e.g. [START_REF] Truscott | Ocean Plankton Populations as Excitable Media[END_REF][START_REF] Hernández-Garcõa | Sustained plankton blooms under open chaotic flows[END_REF][START_REF] Malchow | Spatiotemporal patterns in an excitable plankton system with lysogenic viral infection[END_REF]. It is important to add as well that theoretical studies predict enhancing of ecosystems' stability in case the functional response of predator is of sigmoid type (Oaten and Murdoch, 1975a;Scheffer and De Boer, 1996).

What is the most frequently observed type of functional response for herbivorous zooplankton? The conventional way of determining the type of response consists in laboratory experiments on zooplankton feeding, which are usually carried out in small-sized containers or bottles. A huge amount of literature exists on this topic. Experiments reveal non-sigmoid types (i.e. Holling types I and II) for most herbivorous zooplankton species (e.g. [START_REF] Demott | Feeding selectivities and relative ingestion rates of Daphnia and Bosmina[END_REF]Hansen et al., 1990, Hirst andBunker, 2003;[START_REF] Jeschke | Consumer-food systems: why type I functional responses are exclusive to filter feeders[END_REF]. As a result, a rather strong opinion in the literature is that implementation of Holling III type in plankton models is biologically meaningless [START_REF] Demott | Feeding selectivities and relative ingestion rates of Daphnia and Bosmina[END_REF]Scheffer and De Boer 1996;[START_REF] Murdoch | Plankton abundance and dynamics across nutrient levels: tests of hypotheses[END_REF]. On the other hand, there are a number of papers claiming exactly the opposite, i.e., that this is Holling type III that we should use in modelling of grazing by herbivorous [START_REF] Frost | A threshold feeding behavior in Calanus paczficus[END_REF]Chow-Fraser and Srules, 1992;[START_REF] Gismervik | Prey switching by Acartia clausi: experimental evidence and implications of intraguild predation assessed by a model[END_REF][START_REF] Sarnelle | Type III functional response in Daphnia[END_REF]. More surprisingly, the latter statement is sometimes based on revisiting previous laboratory experiments which had demonstrated non-sigmoid response for a given species (e.g. [START_REF] Sarnelle | Type III functional response in Daphnia[END_REF]. Thus, the question of adequacy of implementation of Holling type III response in models is far from being solved.

I would say, however, that the above discussion on the existence (or non-existence) of type III functional response in feeding experiments might be of little practical importance when one The point is that in all models without explicit vertical space, the grazing rate of zooplankton in the column (per unit volume) is computed as:

( )Z P f E = , ( 1
)
where P and Z are the average (over the column) densities of phytoplankton and zooplankton respectively, f (P) is the functional response found previously in laboratory experiments.

However, the mean-field approach based on (1) provides satisfactory results only in case the vertical distribution of plankton is close to homogenous. On the other hand, it is well known that real vertical distribution of plankton is often highly heterogeneous, especially during phytoplankton blooms and this can seriously affect the behavior of plankton models (e.g. [START_REF] Poggiale | Enrichment paradox induced by spatial heterogeneity in a phytoplankton -zooplankton system[END_REF]. Moreover, herbivorous zooplankton is capable of adjusting its vertical location and feeding at depths of high food concentrations [START_REF] Townsend | Subsurface chlorophyll maxima and vertical distribution of zooplankton in the Gulf of Maine[END_REF][START_REF] Napp | Vertical distribution of marine particles and grazers. II. Relation of grazer distribution to food quality and quantity[END_REF][START_REF] Castro | Vertical distribution of copepods and the utilization of the cholrophyll a-rich layer within Concepcion Bay, Chile[END_REF][START_REF] Seuront | Microscale complexity in the ocean: turbulence, intermittency and plankton life[END_REF]. As such, the grazing rate E in all models with implicit vertical space should be given by E=F/H, where F is the consumption of food in the entire water column, H is the depth of the column. The grazing rate (1) can be re-written in the following way:

( ) ( ) [ ]Z h Z h P g Z Z H F E , = = , (2)
where the ratio ( )

Z H F /
is replaced with the function g [P(h), Z(h)], which is a 'true' functional response, i.e. the amount of food consumed in the column per zooplankton biomass.

In this paper, I show that the functional response g (which I call further the overall functional response) might be of sigmoid type even though the local grazing is described by a nonsigmoid response. Note that the overall response g is not a function in a strict mathematical sense but it is a function of functions (i.e., a functional) since it depends of vertical distributions of phyto-and zooplankton P(h) and Z(h) which are functions themselves. However, in this paper, I will consider the case when ) (P g g =

, where P is the space average density of phytoplankton in

Evidence for Holling type III zooplankton functional response in situ

First, a rather detailed example of emergence of the overall response of sigmoid type from the local non-sigmoid one can be found in Morozov et al. (2008). The ecosystem under study was (Bautista and Harris, 1992). The results are represented in Fig. 1. Fig. 1A shows the overall functional response g for Calanus helgolandicus (large-sized copepod species) plotted versus the average chlorophyll density P in the water column (H=50m). In other words, the specific consumption of food of the whole population of Calanus helgolandicus is plotted versus the average food density in the column. Note that the constructed functional response corresponds to the pre-bloom and the beginning of the bloom periods (March-April 1989). Such a restriction is needed due to the fact that a succession among different algal groups occurs at the peak of the bloom as well as during the post-bloom period (Bautista and Harris, 1992). As such, different algal groups have different nutrition properties and constructing the functional response as a function of total amount of chlorophyll does not make any sense in this case.

To demonstrate that the plotted in Fig. 1A functional response is closer to Holling type III, the data were fitted with the power functions γ α P g = (with Ȗ ≥ 1). The latter fitting can be considered as the first term of expansion of the sigmoid response given by ) 1 /( γ γ β α

P P + ,
which is well adapted in the literature (see [START_REF] Gentleman | Functional responses for zooplankton feeding on multiple resources: a review of assumptions and biological dynamics[END_REF] and the references therein). In section 3.3, some theoretical background will be given for the approximation of the overall response by a power function. The non-linear regression analysis based on the least square method (LSM) has been applied to provide the best parameters (which are given in legends of Fig. 1). The analysis shows that the non-linear responses given by γ α P g = fit the data better than the linear one (for the linear response R 2 =0.52, Ȗ=1; for the nonlinear responses: R 2 = 0.706, Ȗ =2; R 2 =0.707, Ȗ =3). In each case, the parameter Ȗ was fixed at the indicated values. Note that fitting both parameters Į, Ȗ simultaneously results in large values of confident intervals for Ȗ.

As a result, the overall response of the zooplankton population shown in Fig. 1A [START_REF] Schnack | On the feeding of copepods on Thalassiosira partheneia from the Northwest African upwelling area[END_REF]). The possible explanation of such alteration between the responses consists in changing in feeding strategies of zooplankton with an increase of total amount of food in the system (Bautista and Harris, 1992). In particular, the major grazing takes place at depths with high food abundance.

To make the conclusion about the Holling type III response in Fig. 1A more convincing, I plotted as well the clearance rate of zooplankton which is determined by

P g q / =
. Note that the clearance rate of a sigmoid response increases with an increase of food density at low food densities whereas the clearance rate of non-sigmoid response either always decreases (Holling type II) or remains constant for Holling type I (Chow-Fraser and Srules, 1992). Fig. 1B demonstrates an increase of the overall clearance rate q with an increase of food which is in support of sigmoid shape of functional response. Based on the previous approximation of g by power functions, I considered fitting of q as 1 -= γ α P q . This gives R 2 =0 for the linear functional response (Ȗ=1); R 2 = 0.47 for parabolic functional response (Ȗ =2) and R 2 =0.40 for a response parameterized by a cubic parabola (Ȗ =3). Note also that implementation of Akaike information criteria tells us that the linear response has substantially less support compared to responses with Ȗ >1. This takes place for approximation of response function g and the clearance rate q.

Surprisingly, there is a notable scarcity of data in the plankton literature which could allow to construct overall functional responses of zooplankton in situ similar to the one constructed above. This is especially true for data on zooplankton feeding during pre-bloom and the beginning of bloom periods. Note that these are exactly those periods where one can distinguish between sigmoid and non-sigmoid responses since all responses behave in a similar way at high phytoplankton densities, exhibiting saturation. However, the mentioned difficulty can be partially settled in case one is interested in demonstrating only the existence of Holling type III response without plotting the response itself.

In this paper, I suggest a simple straightforward test to demonstrate the emergence of Holling type III response in real plankton communities. To implement the test, one needs to plot

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

8 the ratio of the overall functional responses g 2 /g 1 versus the ratio P 2 /P 1 of the total amounts of food at any two stations. Note that the environmental conditions (e.g. temperature, hydrological regimes, etc) at the two stations should be very close to each other. The results can be easily interpreted with the help of (g 2 /g 1 , P 2 /P 1 ) diagram (see Fig. 2A). In such representation, Holling type I would correspond to the points lying on the dotted line g 2 /g 1 =P 2 /P 1 , i.e. an increase of consumption rate would be proportional to the increase of the amount of food. Points related to

Holling type II response should get inside the domain situated below the dotted line. Note that, the same concerns Holling type III response after saturation. Finally, the domain, corresponding to

Holling type III (before the saturation) is situated above the dotted line since an increase in food amount by the factor of m, would result in an increase of average specific grazing rate by the factor of n>m. Note that we consider that P 2 /P 1 >1, otherwise, the location of domains for Holling type II and III will be the opposite.

Implementation of the above test demonstrates the evidence for Holling type III overall response in different ecosystems (see Fig. 2B). I used literature field data on grazing of copepods in the ocean, making sure that there was no pronounced succession among different algal groups with different nutrition characteristics during the considered periods. One can see from Fig. 2B that for the considered ecosystems, the points in the diagram are situated well above the dotted line,

showing an increase in grazing rate (per unit zooplankton biomass) more rapid than linear. The details regarding the data sets used in Fig. 2B are given in the figure legend. Thus, the responses g(P) of the whole zooplankton populations in the column as functions the total amount of food were of sigmoid type. On the other hand, it is to be emphasized that grazing of each of those species in laboratories shows a non-sigmoid 'local' response.

Modeling the emergence of Holling type III zooplankton functional response

General idea of the method

The central idea of the method implemented below is the following. A plankton model with explicit vertical resolution is used to obtain the overall functional response. Based on the

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

vertical distribution of species given by the model's equations, the specific intake rate of the whole zooplankton population is plotted versus the average density of food in the column. The overall functional responses obtained in such way might be used when constructing plankton models with implicit vertical space.

To compute the intake rate in a model with vertical resolution, one needs to know the distribution of food in the column as well as that of grazers. It is important to emphasize that one needs the vertical distribution of actively grazing zooplankton which can be different from the actual distribution of zooplankton in the column. The point is that feeding of zooplankton organisms such as copepods includes periods of active grazing and non-active periods (see [START_REF] Leising | Copepod foraging and predation risk within the surface layer during night-time feeding forays[END_REF] and the references therein). As such, one needs to take into account the grazing impact only of those animals which are currently in the active phase of a feeding cycle and neglect those which are in the non-active phase of feeding (Morozov and Arshkevich, 2010). Let us denote the density of actively grazing zooplankton at a depth h by Z(h) and suppose that consumption of active grazers at a depth h can be described via a local functional response f h (P(h)). Then the overall response g of the whole population (per unit biomass) will be determined by

( ) 0 0 )) ( ( Z H dh h Z h P f g H h ³ = , (3)
where Z 0 is the average zooplankton density (H Z 0 gives the total biomass of zooplankton in the column per unit area), H is the depth of the column. The shape of the overall response g can be revealed by plotting g versus the space average of food density P .

Note that g might not be a unique function of the average food density P , since it depends on distributions P(h) and Z(h) rather than on P . In case P(h), Z(h) are functions of time given by a dynamical model, the 'number' of model's trajectories will be infinite and this would result in an infinite number of functional responses (3). To settle (at least, partially) this problem, a possible solution can be considering the functional response after all transients die out and the system 3) should be understood as time average species densities. Note, however, that such an approach can provide only a few points for a fixed set of parameters (corresponding to the number of attractors). This, definitely, cannot be used for plotting ( ) P g . To obtain a continuous range of P , a system parameter is to be varied. In this paper, the growth rate of phytoplankton is chosen as a control parameter. This would allow to model an ecologically important case of occurrence of algal blooms arising as a response to an increase of water temperature, light intensity, etc (e.g. [START_REF] Truscott | Ocean Plankton Populations as Excitable Media[END_REF]).

Model formulation

I assume the total amount of zooplankton in the column to be constant; however, the vertical distribution of actively grazing zooplankton Z(h) can vary depending on phytoplankton distribution. The influence of the variation in the total amount of zooplankton on the shape of functional response is discussed in 3.3. The dynamics of phytoplankton in the column is described by the following equation ( ) ( ) ( ) ( )

P Z P P f P P h r h P D t P h - + ∂ ∂ = ∂ ∂ , 2 2 , ( 4
)
where D is the coefficient of vertical turbulent diffusion, r is the growth rate of phytoplankton.

Following Morozov and Arshkevich (2010), I shall assume that the distribution of the actively grazing zooplankton follows the ideal free distribution (i.e. it is proportional to the food density in the column)

o Z P P P Z θ = ) ( . (5)
The total number of actively grazing zooplankton in the population given by H o Z θ is assumed to be constant, where ș is a coefficient (0<ș <1).

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

To parameterize phytoplankton growth rate, the logistic function ( ) , where B h is the integral biomass of phytoplankton located between the water surface and depth h; Ȟ is a positive coefficient [START_REF] Herman | Numerical modelling of diel carbon production and zooplankton grazing on the scotian shelf based on observational data[END_REF]. Note that including absorption of light by water into the generic model ( 4) does not influence the main results (see the end of Section 3.3).

K P r / 1- is used,
After having parameterized the functions in (4), we arrive at the following equation

P Z P P K P P dx P r h P D t P o h β α ν + - 1 • © § - 1 • © § - + ∂ ∂ = ∂ ∂ ³ 1 1 1 exp 2 0 0 2 2 , ( 6
)
where the local functional response is non-sigmoid (Holling type I for ȕ=0 and Holling type II for ȕ0). The coefficient ș is included in Į; the parameters Į, ȕ are considered to be depthindependent; P is the space average density of phytoplankton.

The unit of plankton density is chosen to be ȝg C l -1 . The maximum phytoplankton growth rate approximately belongs to the range 0.1< 0 r < 4 1/day [START_REF] Edwards | Zooplankton mortality and the dynamical behavior of plankton population models[END_REF]. The system's carrying capacity K is rather hard to be estimated. I shall consider K=100 ȝg C l -1 which is consistent with some literature (e.g. [START_REF] Franks | Phytoplankton blooms in a fluctuating environment: the roles of plankton response time scales and grazing[END_REF]. The coefficients describing the local functional response of zooplankton can be estimated as 0.01<α <0.3 1/ȝg C l -1 /day; 0.005< β <0.2 1/ȝg C l -1 [START_REF] Hansen | On the trophic fate of Phaeocystis pouchetii (Hariot). III. Functional responses in grazing demonstrated on juvenile stages of Calanus finmarchicus (Copepoda) fed diatoms and Phaeocystis[END_REF] ). The attenuation constant Ȟ is given by 005 . 0 = ν 1/(m ȝg C l -1 ) [START_REF] Herman | Numerical modelling of diel carbon production and zooplankton grazing on the scotian shelf based on observational data[END_REF]Beckmann and Hense, 2007). I consider the diffusion coefficient to be D= 1 m 2 /d which is in agreement with the data (Beckmann and Hense, 2007). The density of zooplankton in the column (space average density) is assumed to vary between 1< Z 0 < 20 ȝg C l -1 . gives the half-saturation density of feeding). The dashed curves (with the same numbers) represent the corresponding local responses which would take place in case of a well-mixed system. Fig. 3A gives functional responses for large half-saturation densities (including Hoilling type I local response having 1/ȕ =). Fig. 3B shows functional responses for smaller values of half-saturation density.

One can see from the Fig. 3A that for small P , all the overall responses are close to linear. This is because the distribution of food in the column is close to homogeneous (the self-shading is negligible) and the influence of spatial heterogeneity is week. Further, with an increase of P , the overall responses become of Holling type III (concave upwards). This is despite the fact that local functional responses are concave downwards. Increase in ȕ shows (see Fig. 3B) that the range of P corresponding to a concave upwards overall functional response is shrinking. In particular, for a supercritical ȕ, there is no alteration of type between the local and the overall response.

The emergence of Holling III type overall response can be seen as well based on the analysis of the clearance rates ( ) P q , obtained by dividing the functional response by the food density. Murdoch, 1975a;[START_REF] Bazykin | Nonlinear Dynamics of Interacting Populations[END_REF].

To understand the mechanism of emergence of Holling type III response in the model, it is convenient to follow the changes in phytoplankton vertical distribution as the total algal biomass increases. Fig. 4 shows profiles of the absolute and the relative algal densities (Fig. 4A and Fig. 4B, respectively) for different P . One can see that the degree of heterogeneity in algal vertical distribution becomes more pronounced while the total amount of phytoplankton increases. Such spatial structuring of phytoplankton in the model occurs due to the process of self-shading. The distribution of active grazers follows that of food (according to (5),) and the major grazing becomes shifted to depths with higher food density. As a result, the overall grazing rate in the column becomes larger compared to the situation when same amount of food is distributed homogeneously.

In appendix A, the shape of the overall functional response is computed for K>>1 and ȕ<<1. Considering only two first non-zero terms of (A8) gives ( )

( ) 1 • © § + ≈ 12 1 2 3 H P P P g ν α (7)
and for the overall clearance rate ( )

( ) 1 • © § + ≈ 12 1 1 2 2 H P P q ν α . ( 8
)
The stationary vertical profile of phytoplankton is given by (A4). From ( 7) one can conclude that the increase of the overall functional response for intermediate food densities can be approximately described by a power function. This gives some theoretical background for the implementation of non-linear fitting of the real data in section 2 in the form of The difference between the local and the overall functional responses becomes less pronounced in plankton ecosystems in shallow waters. Fig. 5 shows the overall functional responses (the upper row) in a plankton ecosystem with a small depth (H=10m). The corresponding clearance rates are given in the lower row. All the graphs are constructed for the same parameters Į and ȕ as in Fig. 3. The corresponding local functional responses are shown by dashed lines. One can see that that for the local response of type I (ȕ=0), the emergence of the overall response of type III (sigmoid response) still takes place. However, the range of ȕ, which guarantees alteration of types, shrinks substantially compared to situation in deep waters. For instance, in the considered case, the existence of overall Holling type III becomes impossible already for ȕ>0.005 1/ȝg C l -1 . Note also that the prediction that the convex behaviour of the overall response becomes less pronounced in shallow waters may also be obtained from expression (7).

Interestingly enough, for supercritical values of ȕ, the overall functional response in shallow waters becomes undetermined for small P (see Fig. 5B,D). This is because the system has no stationary states with the space average cr P P < < 0 . In this case, a decrease of the growth rate r 0 from large values to a certain critical value r 0cr would result in a monotonic continuous decrease in P up to cr P . A further decrease of r 0 would result in a discontinuous drop of P to zero. While varying parameter r 0 in the opposite direction, the whole picture will be different. For small values of r 0 , the only attractor will be a zero stationary state. A further increase in r 0 will not change the picture (even if r 0 > r 0cr ) until the zero state becomes unstable, which happens for r 0 >μ. The loss of stability results in a discontinuous sudden increase of P to large values. A further increase in r 0 results in a monotonic increase in P . Since μ > r 0cr , there is a region of bistability in the parameter

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

space with two stable stationary states, one of the states being the zero state. Note that bistability and the hysteresis phenomena in the model occur for small H (shallow waters) and for large ȕ.

Analysis of the influence of the total amount of zooplankton (Z 0 ) on the functional responses in the system shows that in shallow waters such dependence becomes more pronounced compared to the situation in deep waters. In particular, an increase in Z 0 results in an increase of overall functional response. Numerical simulation shows that this is a consequence of a more pronounced vertical gradient in phytoplankton distribution for larger Z 0 . However, the difference between overall responses constructed for 1< Z 0 < 20 does not exceed 5% of the value of g.

Finally, I investigated the model's properties in case the absorption of light by water is taken into account. This can be described via a standard way by multiplying the algal growth rate by ) exp( h λ -

, where Ȝ is the light absorption coefficient [START_REF] Herman | Numerical modelling of diel carbon production and zooplankton grazing on the scotian shelf based on observational data[END_REF]. The light absorption coefficient varies within the following range 05 . 0 01 . 0 < < γ

1/m (Herman and Platt,1983;[START_REF] Vila | Modelling spectral irradiance in freshwater in relation to phytoplankton and solar radiation[END_REF]. The resulting algal growth rate becomes

1 • © § - - - = ³ K P dx P h r r h 1 ) exp( 0 0 ν λ . ( 9
)
Numerical simulation shows that adding of absorption of light by water does not result in qualitative changes in behaviour of the overall functional response. Moreover, it facilitates emergence of Holling type III overall response; the increase in the clearance rate with P for Ȝ>0 is observed within larger parameter ranges compared to Ȝ=0. In particular, Holling type III is found for smaller values of half-saturation density (1/ȕ) compared to the system with Ȝ=0. The corresponding graphs for ( ) P g and ( )

P q for 0 > λ
are not shown here for the sake of brevity.

An approximate expression for ( )

P g
in case the algal growth rate is given by ( 9) is computed in Appendix A (obtained for ȕ <<1, K>>1). Considering the first three non-zero terms of (A14) gives (  ) 24) exp( 1 6 1 1 ) exp( 2) exp( 1P H 10) the corresponding clearance rate q Ȝ can be easily obtained by dividing (10) by P . Note that expression (10) transforms itself into (7) when Ȝĺ0. Comparison of the overall responses ( ) P g and ( ) P g λ reveals that ( ) ( )

1 • © § - + + + - + ≈ 3 2 3 2 2 1 ) exp(
H H P H P H H H P g λ ν λ ν λ ν λ λ α λ , (
P g P g > λ .
Also, the increase in overall clearance rates with P (characterized by ( )

P q λ ′
) becomes larger for the system with Ȝ>0. This can be interpreted as an enhancing of system's stability (Oaten and Murdoch 1975a). The property that ( ) ( )

P g P g > λ can
be explained by the fact that the gradient in the phytoplankton vertical distribution becomes more pronounced when taking into account absorption by water. Since the active grazers follow the distribution of food, the overall consumption of food in the system with Ȝ>0 becomes higher compared to the situation with Ȝ=0.

Discussion

Dynamics of plankton models is sensitive to the choice of type of functional response of herbivores. Implementation of Holling type III (sigmoid) response would enhance stability of models. There is an opinion in the literature (based on feeding experiments in laboratories) that the use of Holling III type in plankton models is biologically meaningless [START_REF] Demott | Feeding selectivities and relative ingestion rates of Daphnia and Bosmina[END_REF]Scheffer and De Boer, 1996;[START_REF] Murdoch | Plankton abundance and dynamics across nutrient levels: tests of hypotheses[END_REF]. In this paper, I challenge this opinion by providing field evidence and theoretical background for the existence of Holling type III zooplankton functional response. I consider the overall functional response of zooplankton in the water column, i.e. the overall intake rate per zooplankton biomass as a function of average food density. Note that this is the overall functional response which is used in all models with implicit vertical space (mean-field models).

My literature search on data on zooplankton feeding in situ has revealed that zooplankton populations in the water column (the overall response) often exhibit a sigmoid behaviour when the local functional responses, observed for the same species in laboratory, are non-sigmoid (i.e. pronounced with an increase of total amount of phytoplankton in the column (due to self-shading of algae). Grazing of zooplankton becomes shifted to layers with higher food abundance and the total intake rate becomes larger as it were the case of the homogeneous distribution of the same amount of food in the column (a well-mixed system). As such, an increase of the total amount of food in the system by the factor of n>1would result in an increase of the intake rate by the factor of m (m>n). Interestingly, a more detailed analysis of causes of emergence of sigmoid response in real ecosystems (found in references in Fig. 2) indicates that the above mechanism might provide a sound explanation of this phenomenon.

The model predicts emergence of Holling type III response to be facilitated in deep waters ecosystems; the light attenuation should be also well pronounced (these conclusions can be drawn directly from analysis of expressions ( 7), (9) for the overall response). This is related to the fact that in shallow waters the spatial distribution of food becomes more homogeneous and the overall functional response is closer to the local one which is non-sigmoid. Another important requirement is that the half-saturation density in the local response should be sufficiently high.

This signifies that the overall response of sigmoid type would arise for large-sized herbivorous,

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

18 usually characterized by smaller values of ȕ. Interestingly, this conclusion is in a good agreement with data on zooplankton feeding in situ. For example, the emergence of Holling type III response shown in Section 2 (see Fig. 1), based on data from (Bautista and Harris, 1992) for large copepod species, is not observed for small-sized species, which exhibit the overall functional response of Holling type II.

The emergence of the overall Holling type III response, shown in the paper, is of vital importance for mathematical modelling of plankton dynamics. In particular, this concerns construction of models without explicit vertical space. The results of the paper suggest that the functional response of the whole zooplankton population in mean-field models might be described by type III response, provided the depth of the column is sufficiently large. This would enhance the stability properties of the system (Oaten and Murdoch, 1975a;Scheffer and De Boer, 1996).

The choice of a concrete parameterization of the overall functional response to be used in meanfield plankton models is a more complicated issue. One can implement expressions ( 7), ( 10) derived analytically, which provide good approximation of ( ) P g or small ȕ. However, in case one is interested mostly in qualitative properties of models, the 'classical' parameterization of Holling type III given by ) 1 /( γ γ β α P P + (Ȗ>1) might be implemented as well by providing a rough approximation of graphs in Fig. 3. In this case, the half-saturation constant (1/ȕ ) 1/Ȗ will have a different meaning from the one measured in laboratory since the parameter ȕ in the overall response would include not only the life traits of species but some characteristics of the habitat as well (e.g. depth, light attenuation, vertical diffusion, etc). On the contrary, in ecosystems with shallow waters the implementation of Holling type II for small P in mean-field models becomes more appropriate. In particular, the absence of functional response in model describing shallow waters (see Fig. 5B,D) shows that zooplankton cannot control phytoplankton growth at low algal density. This is exactly what would predict a mean-field model with type II functional response.

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

In this paper, the overall functional response was considered to be a function of total amount of food in the column, i.e. ) (P g g = . However, as it follows from (2), different vertical distributions P(h) and Z(h) might correspond to the same average biomasses P , Z making it questionable computing the grazing rate in terms of the product Z P g ) ( . However, field observations, especially in nutrient-rich ecosystems, reveal that the consumption rate of zooplankton (per unit biomass) in the column can be approximately considered as a function of the total amount of food (e.g. Morozov et al., 2008, see also Fig. 1 in this paper), i.e. by ) (P g g = .

This can be understood from the fact that the profile P(h) is a result of interplay of different biotic and abiotic factors, thus it cannot be arbitrary. In particular, in eutrophic ecosystems there is a strong feedback between vertical distribution of phytoplankton and its total amount P in the column (via the light attenuation). A typical scenario consists in an increase in the gradient of algal spatial distribution when the total amount of phytoplankton is being increased (e.g. [START_REF] Raymont | Plankton and Productivity in the Oceans[END_REF][START_REF] Ohman | The demographic benefits of diel vertical migration by zooplankton[END_REF]. Such property can be seen from a generic model (6) as well (see Fig. 5). I should say, however, that in oligotrophic ecosystems, there is a possibility of alternative density profiles of phytoplankton (i.e. the existence of surface and deep chlorophyll maxima) for the same environmental characteristics [START_REF] Ryabov | Vertical distribution and composition of phytoplankton under the influence of an upper mixed layer[END_REF]. Probably, dynamics of such bi-stable heterogeneous ecosystems can be described adequately only via spatially explicit modelling approach.

The model's predictions are based on an important assumption about the foraging strategy of predators. The distribution of actively feeding zooplankton (i.e. the organisms which are currently in the active stage of feeding) in the column was considered to be the ideal free distribution. In reality, zooplankton shows more complex food searching behaviour by optimizing its feeding under different constrains (e.g. temperature, presence of high predators, competition among zooplankters of same trophic level, etc). Incorporation of the mentioned constrains would result in vertical profiles of zooplankton deviating from the ideal free distribution law (e.g. Liu et , 2003;[START_REF] Lampert | Vertical distribution of zooplankton: density dependence and evidence for an ideal free distribution with costs[END_REF][START_REF] Leising | Copepod foraging and predation risk within the surface layer during night-time feeding forays[END_REF]. However, this should not be considered as a serious restriction for applications of the obtained modelling results. The point is that the overall sigmoid response in model ( 6) is observed when the ideal free distribution of predators is not really 'ideal'. Simulations show for small perturbations of the ideal free law (5), the increase of P results in a larger than linear increase of functional response, i.e. in a type III response. In particular, this occurs for spatially periodical perturbations having random amplitude and random phase (up to 20-50% of the values of Ș). Another important issue is that including the mentioned constrains, when modelling the optimal foraging strategy of herbivores, would be ecosystems' specific. On the contrary, positive correlations between vertical location of active feeders and food density in the column have been reported in a large number of plankton ecosystems in different geographic regions (i.e. it is a generic property of most plankton ecosystems). Thus, the free ideal law can be considered as zero order approximation of vertical distribution of active feeders for implementation in generic plankton models.

Also, I would like to mention that the assumption about the ideal free distribution of active grazers in the column is not crucial to obtain a sigmoid overall functional response. For instance, a sigmoid-shaped overall functional response has been recently found in a model ecosystem (Morozov and Arashkevich, 2008) in case the rate of active vertical displacement of zooplankton was proportional to the gradient of food. Note that, the resultant vertical distribution of zooplankton in the cited paper was not the ideal free distribution. However, the main mechanism of emergence of Holling type III response was similar to the one found in the current paper.

Finally, it is to be emphasized that the mechanism of emergence of Holling type III zooplankton functional response suggested here is somewhat different from those given earlier in the ecological literature (both in aquatic and terrestrial ecosystems). For instance, the most frequently used explanations of the sigmoid-like shape of the functional response are: (i) the existence of alternative food (Oaten and Murdoch, 1975b;[START_REF] Akre | Switching and sigmoid functional-response curves by damselfly naiads with alternative prey available[END_REF][START_REF] Gentleman | Functional responses for zooplankton feeding on multiple resources: a review of assumptions and biological dynamics[END_REF][START_REF] Elliott | Prey switching in four species of carnivorous stoneflies[END_REF]), (ii) refuge for prey [START_REF] Holling | The functional response of vertebrate pred-ators to prey density[END_REF][START_REF] Taylor | Predation[END_REF][START_REF] Poggiale | Predator-prey models in heterogeneous environment: emergence of functional response[END_REF] [START_REF] Holling | The functional response of predators to prey density and its role in mimicry and population regulation[END_REF][START_REF] Luck | Prey or host searching behaviour that leads to a sigmoid functional response in invertebrate predators and parasitoids[END_REF][START_REF] Van Leeuwen | How population dynamics shape the functional response in a one-predator-two-prey system[END_REF]. None of them is needed to provide a sigmoid response in model ( 6).

Arising of a sigmoid functional response due to the heterogeneity of space is less addressed in the literature. Recently, it was shown that prey patchiness and predator aggregation might result in a sigmoid overall functional response when the local response is non-sigmoid [START_REF] Nachman | A functional response model of a predator population foraging in a patchy habitat[END_REF]. The ecosystem under study was an acarine predator-prey system. In the current paper, a sigmoid-like overall response arises via a close scenario: the ideal free distribution can be considered as aggregation of predator. However, an important difference of plankton ecosystems is the existence of a strong feedback between the spatial distribution of prey (algae) and its total amount in the system (via light attenuation). In this sense, nutrient-rich plankton ecosystems are rather unique among other ecosystems. They combine the light-mediated aggregation of prey distribution (see Fig. 4) and the existence of fast-moving predator capable of regulating grazing load within the whole habitat due to possibility of fast displacement of organisms in the vertical direction. It is quite possible, however, that some other non-plankton ecosystems with similar properties of spatial predator-prey interactions would result in the emergence of an overall Holling type III response according to the reported mechanism. After differentiating both sides of the equation with respect to h we obtain ( ) ( )

h P h P ′ = -2 ν . ( A3
)
Integration of (A3) gives the stationary profile of phytoplankton ( )

C h h P + = ν 1 , ( A4
)
where C is a constant of integration. Consider the space average of both sides of (A4)

1 • © § + = C C H H P ν ν ln 1 . (A5)
The stationary profile of phytoplankton can be rewritten as a function of P ( )

1 • © § - - = H P H h P ν ν exp 1 1 . (A6)
The overall functional response is given by ( )

dh P h P H P g H ³ = 0 2 ) ( 1 α . ( A7
)
By substituting (A6) into (A7) we obtain after simplification ( ) ( ) ( )

( ) P H H P H P P g 2 2 exp exp ν ν ν - - + = . ( A8
)
Taylor expansion of (A8) gives the overall functional response in the following form

( ) ( ) ( ) ( ) 1 • © § + + + + = Κ ! 8 2 ! 6 2 ! 4 2 6 7 4 5 2 3 H P H P H P P P g ν ν ν α . ( A9
)
The clearance rate can be easily computed by dividing (A9) by P ( ) One can see that q is an increasing function of P since its first derivative is positive. Numerical simulation shows that expression (A9) provides a satisfactory approximation of the overall response for a finite values of carrying capacity K (up to K=50) and small values of ȕ>0. This approximation can be used as well in case the vertical diffusion is non-zero D>0 (but small). In particular, it can be used for D=1.

( ) ( ) ( ) 1 • © § + + + + = Κ ! 8 2 ! 6 2 ! 4 2 1 6 6 4 4 2 2 H P H P H P P q ν ν ν α . ( A10
Now, let us consider the situation when the absorption of light by water is taken into account. The growth rate of phytoplankton is given by ( 9). The expression for the overall functional response can be derived in a similar way as above. The differential equation for the vertical profile of phytoplankton is the following ( ) ( )

h P h P ′ = - - 2 ν λ . ( A11
)
Integration of (A11) results in ( )

1 ) exp( ) exp( - = x C x C h P λ ν λ . ( A12
)
where C is a constant of integration. By applying the same technique as before, the vertical profile of phytoplankton can be expressed via the space average algal density P as

( ) ( ) ( ) ( ) ( ) ( ) (
)

H P H h H h P h H H P H h h P ν λ λ ν ν ν λ λ + - - - - - + + + - - - - = ) ( exp ) ( exp ) ( exp 1 1 ) exp( ) ( exp . (A13)
The overall functional response is given by integration of (A13) over the column. However, the analytical expression for ( ) P g λ results in a rather cumbersome analytical formula. The Taylor expansion gives the following expressions for ( )

P g λ ( ) ( ) Κ + + + = 3 3 2 2 1 P a P a P a P g λ , (A14)
where Simulation shows that the three first terms in (A14) give a good approximation of ( ) P g λ in the considered parameter range.

( ) ( ) ( )
1 ) exp( ) exp( 1 ! 2 2 1 2 1 2 1 2 - + + = + + H H H n a n n n λ ν λ λ α , ( )! 1 2 1 2 2 2 + = - n H a n n n ν λ α .

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

1. (A) The overall functional response g of copepods (Calanus helgolandicus) measured in situ.

The data are taken from (Bautista and Harris, 1992). The fitting functions are: (1) linear function 2. (A) Determining the type of functional response based on (g 2 /g 1 , P 2 /P 1 ) diagram, where P 2 /P 1 is the relative increase of total amount of food and g 2 / g 1 is the relative increase in the specific consumption rate of a zooplankton population. ȕ 3 =0.14 1/ȝg C l -1 ). Note that for large ȕ the overall functional response in shallow waters becomes undetermined for small P , for details see the text. The other system's parameters are the same as in Fig. 3.

  the rate of food intake by the whole zooplankton population in the column, as a function of total amount of food, often exhibits a sigmoid behaviour, instead of a nonsigmoid one postulated previously based upon laboratory experiments. This conceptual discrepancy is due to the ability of zooplankton to feed mostly in layers with high algal density. I propose a generic model explaining the observed alteration of type between the overall and the local functional responses. I show that emergence of Holling type III in plankton systems is due to mechanisms different from those well known in the ecological literature (e.g. food search learning, existence of alternative food, refuge for prey).

  zooplankton feeding in real ecosystems. This is especially true for plankton models without vertical spatial resolution, dealing with average over the column species densities.

  The correctness of such assumption will be discussed in section 4. Note that such assumption is made in all models without explicit spatial resolution. As a result, the grazing rate of zooplankton is given by e. as in expression (1). However, the laboratory functional response (or the local functional response) ( )oftype of response, in particular, it might result in a Holling type III overall response. As a part of the current study, I have conducted a rather comprehensive literature search on patterns of zooplankton feeding in situ. The results show the evidence of alteration of type of the response between the overall ( ) P g and the local ( ) P f functional responses. In particular, by applying a straightforward test, I demonstrate the emergence of Holling type III overall response for herbivorous zooplankton for a number of ecosystems. Surprisingly enough, sigmoid overall responses are obtained for species usually exhibiting non-sigmoid responses (local responses) in laboratories. Such discrepancy is not due to some errors in previous laboratory studies, but is a result of complex foraging behavior of zooplankton in real ecosystems. Further, I suggested a simple generic plankton model explaining the observed alteration of type of functional response. The model shows emergence of Holling III type for the overall response ( ) P g with a non-sigmoid local response of zooplankton f (P) within a large range of parameters. The obtained Holling type III becomes possible due to a mechanism different from those which are well known in the literature to be responsible for a sigmoid-shaped response (e.g., learning in food searching, existence of alternative food, refuge of prey). The model predicts the crucial role of light attenuation by phytoplankton and the depth of the column in the shaping of the zooplankton functional response.

  Sea. Here I provide another example of such alteration of types based on field data from

  the attractors. In case the attractors are not stationary states but non-damped oscillations (regular or chaotic), the values of P(h), Z(h) in (

  with K being the carrying capacity. Thus, depletion of nutrients in the column is not modelled explicitly. The vertical heterogeneity in phytoplankton growth rate in the model arises due to light shading. I considered the scenario when the attenuation of light in the column is due to selfshading of phytoplankton in deep layers by phytoplankton of surface layers. This can be parameterized by multiplying the algal growth rate by ) exp(

  Results of modellingAnalysis of model (6) shows that the shape of the overall zooplankton functional response largely depends on the column's depths H. Thus, I shall consider here separately functional responses in deep and shallow waters.

Fig. 3

 3 Fig.3shows the overall functional responses (the upper row) in a plankton ecosystem with

  Fig.3C,D show the clearance rates constructed based on responses in Fig.3A,Brespectively. For small ȕ the overall clearance rates exhibit intervals of increase, i.e. ( ) P q′ . An increase of clearance rate with P can be considered as one of criteria of the emergence Holling type III response. Note that theoretical studies predict enhancing stability of predator -prey the clearance rate is an increasing function of the food density(Oaten and

  in Fig.3have been obtained by varying the algal growth rate r 0 for fixed values of Z 0 . Different value of total amount of zooplankton Z 0 changes the overall responses. Interestingly enough, such changes turn out to be rather small within the considered less than 1-2% of the value of g). This is mostly due to the fact that the vertical diffusion (D = 1 m 2 /d) does not play an important role in shaping the vertical profiles of species in the model. Note that in model (6) without diffusion, the overall functional response does not depend on the amount of zooplankton in the system.

  or II). Such behaviour has been found in different ecosystems (see Fig.2). The observed alteration of types between the overall and the local functional responses was a result of interplay between heterogeneous food distribution in the column and active zooplankton foraging behaviour. As such, it would be preferable to implement type III response in the mean-field plankton models even though non-sigmoid types have been unambiguously reported in laboratories.I proposed a simple generic model to explain the observed emergence of a sigmoid overall response. The model shows a surprising discrepancy between the behaviour of local response and that of the overall response. In particular, a Holling type III overall response is observed in case the local functional response is non-sigmoid. The model predicts the following generic mechanism of emergence of a sigmoid response. The gradient of vertical algal distribution becomes more

  2009), (iii) learning in food searching by predator

(

  B) Evidence for the sigmoid functional response in real plankton communities. The points are situated well above the dotted line g 2 / g 1 = P 2 /P 1 which supports a sigmoid-like behavior of the responses. The data represent feeding of marine copepods in the following ecosystems: (1,2) Barents Sea,May and July 2003[START_REF] Lampert | Vertical distribution of zooplankton: density dependence and evidence for an ideal free distribution with costs[END_REF] (Morozov et al, 2008); (3) coastal waters off Plymouth,March-April 1989 (Bautista and Harris, 1992); (4) coastal waters off Plymouth, March-April 1990(Bautista et al, 1992); (5) Faroe Shelf, May 2004[START_REF] Debes | Seasonal variability in copepod ingestion and egg production on the Faroe shelf[END_REF]; (6) Western Mediterranean, June 1986[START_REF] Saiz | Pigment gut contents of copepods and deep phytoplankton maximum in the Western Mediterranean[END_REF]); (7) Bering Sea, May 1979[START_REF] Dagg | The feeding, respiration, and excretion of zooplankton in the Bering Sea during a spring bloom[END_REF]; (8) Eastern Barents Sea (Arashkevich, unpublished data) The bars in the figure show error-bars.3. (A,B)The overall zooplankton functional responses in an ecosystem with a pronounced depth (H=100m). The overall responses are represented by bold curves; the dashed curves (with the same numbers) show the corresponding local responses which would take place for homogeneous vertical distribution of plankton (a well-mixed system). Different curves are constructed for different ȕ. Fig.3Ashows functional responses for large half-saturation densities in the local 1 =0; ȕ 2 =0.02; ȕ 3 =0.04 1/ȝg C l -1 ); Fig.3Brepresents functional responses for small half-saturation densities (ȕ 1 =0.08; ȕ 2 =0.1; ȕ 3 =0.14 1/ȝg C l -1 ). The other system parameters are: K=100 ȝg C l -1 ; Į=0.3 1/ȝg C l -1 m ȝg C l -1 ); D= 1 m 2 /d; Z 0 = 5 ȝg C l -1 . (C,D)The overall zooplankton clearance rates (bold curves) obtained for the constructed above overall functional responses. The dashed curves (with the same numbers) show the local clearance rates. 4. (A) Variation in vertical profile of phytoplankton in the model with an increase in the total amount of phytoplankton P . (B) The corresponding variation in the relative vertical density of phytoplankton P P / . The curves in (A,B) are constructed for the following values of P : -1 ; ȕ=0.02 1/ȝg C l -1 . The other system parameters are the same as in Fig.3.5. The overall zooplankton functionalresponses (A,B) and the clearance rates (C,D) in a shallow ecosystem (H=10m). The overall responses (or the overall clearance rates) are represented by bold curves; the dashed curves (with the same numbers) show the corresponding local responses and local clearance rates. Different curves are constructed for different ȕ. Fig.3A shows functional responses for large half-saturation densities in the local response (ȕ 1 =0; ȕ 2 =0.02; ȕ 3 =0.04 1/ȝg C l -1 ); Fig.3B represents functional responses for smaller half-saturation densities (ȕ 1 =0.08; ȕ 2 =0.1;

  Figure1

Acknowledgements

I highly appreciate Dr. S.V. Petrovskii (University of Leicester) and E.G. Arashkevich (Shirshov Institute of Oceanology, Moscow, Russia) for a careful reading and comments. Also, I am very grateful to both anonymous referees for nice words regarding my paper and their suggestions to improve the manuscript.

Appendix A

Here we derive analytically an expression for the overall functional response in case when K<<1 and ȕ<<1, i.e. when the local functional response is linear or close to linear.

Let us first consider the situation when attenuation of light is due to self-shading of algae.

The model equation is given by (6). To avoid cumbersome formulae, we shall consider the case when the diffusion term in ( 6) is much smaller compared to the growth rate and the grazing rate.

Neglecting the diffusion term and setting