domain: shs sdv

text:
Back to the Root: The Coming of Age of Cementochronology

Stephan Naji and William Rendu Tremendous progress has been made since the 1950s when the first applications of dental cementum analysis were implemented in faunal and human remains (Chapter 1). Today, with a better understanding of cementum biology, robust, replicable protocols support promising innovative research in cementochronology. This chapter summarizes the key advances presented in this volume. We will also address recurring issues that could confuse those who are not familiar with the cementum literature. Our hope for the future is to continue expanding our understanding of cementum from a biological perspective, optimizing the various protocols, and testing broad anthropological hypotheses from the fossil record to forensic contexts.

Part I: The Biology of Cementum

Many aspects of cementum biology are confusing to nonbiologists interested in using acellular cementum's circannual growth pattern in cementochronology (Introduction, Chapter 1). The first part of the volume addressed these issues, and we will summarize what is consensual today among specialists.

Biology

Two main types of cementum are present on teeth roots. Acellular cementum (AEFC) covers cervical root surfaces, and cellular cementum (CIFC) covers apical and furcation regions, even though variations exist between species in the proportion's distribution of AEFC or CIFC along the root. For example, only CIFC is present in horses [START_REF] Sahara | Development of coronal cementum in hypsodont horse cheek teeth: Coronal cementogenesis in horse cheek teeth[END_REF]. Human disease and experimental animal models are now providing new insights into these two cementum types' functions as growth and maintenance processes, respectively (Chapters 2-5). To further our knowledge of cementum, we argue (Introduction) that cementogenesis should be considered within a broader evolutionary framework describing the selection for the thecodont dental attachment morphology probably as early as 200 million years ago [START_REF] Newham | Reptile-like physiology in Early  Jurassic stem-mammals[END_REF]. Empirical evidence suggests that AEFC is under tight genetic control, with a constant growth rate within individuals and a consistent circannual deposit pattern (Introduction,. Cellular cementum provides the complementary adaptive system for tooth attachment against environmental influences, primarily biomechanical loadings and pathological insults (Chapters 3,5). Regarding AEFC composition and structure, clues to what we see when observing a semi-thin undemineralized section under transmitted light microscopy (TLM) start to emerge in Chapters 6, 7, and 14, particularly in terms of growth mineralization variation. The recent development of virtual 3D analysis combining x-ray fluorescence and x-ray diffraction mapping with synchrotron radiation microbeams was precise enough to follow individual increment oscillations across the cementum width. Independent experiments on several species of ungulates, primates, and humans evidenced the regular oscillations of calcium (Ca), phosphorous (P), zinc (Zn), strontium (Sr), and hydroxyapatite crystals (cAp) with seasons and revealed annual cementum band structure. Microtextural variation is also present across individual bands, demonstrating that the overall collagen fibril orientation undergoes subtle changes with the season. The "growth zones" creating wider and lighter (more translucent) increments should probably be interpreted as higher mineralization periods (higher mineral to organic ratio) but loosely organized. In contrast, the dark (opaque) "arrest lines" correspond to more densely packed but less mineralized events. A more comprehensive test of this hypothesis is already underway by some of the volume contributors.

Cementum: A Biological Archive of Life-History Events The most promising research area in cementochronology is probably the potential to identify life-history events precisely. The recognition that cementum growth variations can be potentially affected by a variety of "stressors" such as pathologies (environment), parturitions (development), or fractures (traumatic) has been empirically observed in many mammal species by several authors in the past twenty years. This variation opens the door to exciting life-history studies (Chapters 7, 8).

The consensus is that at least four signals -ontogeny, phylogeny, biomechanics, and environmentare reflected in the micromorphology of skeletal and dental tissues deposited at any time during the life of an individual [START_REF] Padian | Research applications and integration[END_REF][START_REF] Temple | Bioarchaeological evidence for adaptive plasticity and constraint: Exploring lifehistory trade-offs in the human past[END_REF].

The environment part refers to the direct effects on growth and tissue expression caused by periodic, cyclical (seasonal variations), and unique environmental events or stress in a single individual (Jylhävä et al. 2017). To address this topic reliably, broader questions must be tackled first, such as, What variables affect individual cementum deposits' width or how is individual frailty expressed in cementum growth variation throughout life? Today, we can only identify some of these changes in acellular cementum growth pattern in documented samples (Chapter 14). In other mammals, especially bears, these observations are routinely done [START_REF] Matson | A method for dating tetracycline biomarkers in black bear cementum[END_REF]. More valid protocols will be needed to translate these initial observations on documented samples to archaeological samples.

Other Considerations on Cementum Variations

We would like to address some confusing points that are lingering in the anthropology literature. The dietary hypothesis. The change in food coarseness, and consequently in biomechanical loading to process the seasonally soft (summer) and tough (winter) food items in wild populations, and its influence on AEFC composition and optical properties specifically [START_REF] Lieberman | The biological basis for seasonal increments in dental cementum and their application to archaeological research[END_REF], does not seem to be a robust explanatory mechanism across the large numbers of mammals that have been tested, including humans. Empirically, we know of several contexts that directly disprove this hypothesis: (1) Impacted teeth show a comparable amount of AEFC deposits than erupted teeth [START_REF] Azaz | Correlation between age and thickness of cementum in impacted teeth[END_REF][START_REF] Azaz | Aging of tissues of the roots of nonfunctional human teeth (impacted canines)[END_REF][START_REF] Nitzan | The effect of aging on tooth morphology: A study on impacted teeth[END_REF][START_REF] Rai | Effect of nutrition on coronal displacement of cementum in impacted teeth[END_REF] ; (2) vestigial red deer canines also exhibit the same contrasted growth pattern, even though they are not functional for nutrition (unpublished data); and (3) hibernating species display similar seasonal banding even though they do not consume food during winter months [START_REF] Grue | Review of the formation of incremental lines in tooth cementum of terrestrial mammals[END_REF][START_REF] Helm | Circannual rhythms anticipate the Earth's annual periodicity[END_REF][START_REF] Klevezal | Recording Structures of Mammals: Determination of Age and Reconstruction of Life History[END_REF]. Also, genetic analyses suggest that masticatory forces are probably not a driving factor (Chapter 4). Finally, this food coarseness hypothesis would predict that significant changes in the diet, such as the advent of farming (more reliable and predictable softer food) and the industrial revolution (quasi-complete homogenization of diet with processed food), should result in a clear decrease in cementum contrast and growth worldwide. This is not the case [START_REF] Klevezal | Recording Structures of Mammals: Determination of Age and Reconstruction of Life History[END_REF][START_REF] Naji | Cementochronology, to cut or not to cut?[END_REF]. Similarly, the same prediction should apply to wild versus domesticated animals. However, no noticeable changes in AEFC contrasts are visible in animal studies or across human populations [START_REF] Klevezal | Recording Structures of Mammals: Determination of Age and Reconstruction of Life History[END_REF]. If the diet influences the cementogenesis of AEFC, a more likely hypothesis would be that the nutritional status (i.e., the availability of nutrients related to seasonal food) influences growth and mineralization processes, mainly via dietary metabolic pathways. Sporadic deficiencies will be accounted for through energy allocation away from growth processes and toward maintenance, probably throughout the skeleton [START_REF] Bogin | Chapter 11 -The evolution of human growth[END_REF][START_REF] Temple | Bioarchaeological evidence for adaptive plasticity and constraint: Exploring lifehistory trade-offs in the human past[END_REF]. This is essentially the hypothesis we are testing when trying to identify "stress events" in cementum layers. Overall, outside of chronic conditions or genetic defects, nutrition is not expected to affect AEFC deposition in terms of the number of layers (probably under stricter genetic control) but potentially in terms of the mineralization rate or density. However, we would argue that the dietary hypothesis may be valid regarding CIFC growth. We hypothesize that accurate CIFC age-at-death estimates in specific taxa (Chapter 7; [START_REF] Matson | A Laboratory Manual for Cementum Age Determination of Alaska Brown Bear PM1 Teeth[END_REF]) convincingly illustrate its regular and continuous adaptive function under specific conditions. When the feeding pattern and consequent biomechanical pressure are regular, with strong cyclic changes in food coarseness and nutritional value, CIFC growth can mimic the AEFC circannual deposition and be used for age estimation. The seasonal adaptation hypothesis. The dietary hypothesis might be a covariant consequence of the dominating genetic control synchronized with a circannual periodicity (Introduction). This hypothesis is supported empirically for several species the northernmost latitude exhibiting the most contrast in growth layers and, inversely, for the tropical/equatorial samples showing the least contrast (Klevezal' 1996, 101). For example, in tropical regions, cementum layers are more annual than seasonal and are more complex to identify than in temperate areas. Likewise, growth layer formation is inverse in the two hemispheres (Klevezal' 1996, 91) and can shift for a displaced population between the two hemispheres. Also, hibernating species still produce seasonal increments [START_REF] Klevezal | Recording Structures of Mammals: Determination of Age and Reconstruction of Life History[END_REF]Klevezal' & Kleinenberg 1967). Overall, the complexity and variability of cementum growth layers are tied to seasonality gradient (i.e., the more seasonal, the more regular and predictable growth deposits can be identified). Today, the metabolic and endocrine pathways are understood at a much finer level. Some studies have elegantly linked bone growth to nutritional status and metabolic fitness on several ruminant species (Kohler et al. 2012). In summary, during the unfavorable seasons (dry/winter), ruminants are subject to cyclical variation in body core temperature, resting metabolic rate, and associated changes in hormonal levels (especially plasma insulin-like growth factor-1 or IGF-1). The consequence of this energy-conserving strategy is the formation of arrest growth lines. During the good season (humid/summer), the opposite happens (e.g., peak metabolic and IGF-1 levels, enhanced thyroid activity), which stimulates bouts of intense skeletal growth (Kohler et al. 2012). These cyclic growth patterns are synchronized with food availability and nutritional profile, which will dictate, in part, fat storage. Growth arrest and peak are, therefore, concurrent with hormonal and metabolic nadirs and zeniths, respectively. From this model, we predict that the number of circannual AEFC increments is tightly controlled genetically. However, the quality (i.e., increment width, mineralization density, collagen matrix organization) may vary according to external factors or stressors.

Part II: Protocols

The main hurdle in cementochronology and one that we have argued is probably the underlying reason the method never became more ubiquitous [START_REF] Naji | Cementochronology, to cut or not to cut?[END_REF] is the necessity to implement a histology protocol. Laboratory access and a cementumspecific protocol are required to optimize the observations of AEFC increments. Furthermore, numerous researchers have emphasized the critical importance of training. As Buikstra discusses in this volume (Chapter 1), standardization of a protocol has been overdue. This is no longer the case. For freshly extracted teeth with high collagen content, observing decalcified paraffin-embedded, stained ultrathin sections (anywhere between 5 to 15 μm depending on the species) is the most accurate option today [START_REF] Matson | A Laboratory Manual for Cementum Age Determination of Alaska Brown Bear PM1 Teeth[END_REF]. Because ultrathin sections are a standard histology procedure in a forensic context, decalcified thin sections have also been used in humans with high success (Chapter 1: Table 1.1). However, for dried teeth (e.g., archaeological, osteology collections, fossils), this process will destroy the entire root. Undemineralized, embedded semi-thin ground sections (around 100 μm) are the only viable option in archaeological contexts, except for virtual 3D reconstruction.

Optimizing Thin-Section Production

For humans, Wittwer-Backofen and colleagues [START_REF] Wittwer-Backofen | Age estimation using tooth cementum annulation[END_REF][START_REF] Wittwer-Backofen | Tooth cementum annulation for age estimation: Results from a large known-age validation study[END_REF][START_REF] Wittwer-Backofen | Age estimation by tooth cementum annulation: Perspective of a new validation study[END_REF] pioneered the standardization process, upon which Chapters 9-12 have improved. None exists for animals. These chapters propose systematic discussions on thin-section preparation variables and imaging parameters (Chapters 11,12). Consequently, Chapters 1, 9, and 11 propose speciesspecific protocols for human and ungulate teeth. However, though the methods deliver consistent results for experienced specialists, these studies also clearly underline the need for extensive training.

Season-at-Death Validation Studies

In wildlife populations and zooarchaeology samples, estimating the season of death has been implemented reliably for the past forty years (Chapter 1). In humans, seasonat-death estimation is only starting to be explored successfully (Chapter 13;[START_REF] Wedel | Determination of season at death using dental cementum increment analysis* †[END_REF]. In cementochronology, one of the central differences between humans and other mammals is that humans are long lived and thus accrue many cementum bands in a nonlinear pattern [START_REF] Jylhava | Biological age predictors[END_REF]. Consequently, in addition to the relatively small size of human teeth, the last deposited cementum band is challenging to identify optically. [START_REF] Wedel | Determination of season at death using dental cementum increment analysis* †[END_REF] was the first to develop a method in a documented collection of freshly extracted teeth for forensic season-at-death purposes and then applied it to archaeological remains [START_REF] Wedel | Using dental cementum increment analysis to estimate age and season of death in African Americans from an historical cemetery in Missouri[END_REF]. Chapter 12 expands cementochronology to subadults teeth, while Chapter 13 tests with success a comparable protocol on US army personnel.

However, there is a statistical sampling issue with season-at-death estimations that distinguished animal and human studies. In animal analyses, when the last band is a growth zone (summer), its large width allows for two to four subdivisions depending on the author or species (Chapters 17 and 18). Statistically, four outcomes are a critical difference over a binary outcome (summer/winter), which is the only option in humans with a thin, last band. In this scenario, the estimate would require a large sample to support its validity statistically. So far, human methods rely only on a binary (light/ dark) outcome. More tests are needed to exclude random results alone and validate the prediction of season at death.

Cementum Growth Variations and Stress Markers

Several researchers are trying to systematize observations from reference collections on animal and human remains to identify the impact of "stressors" on cementum growth reliably [START_REF] Edinborough | Mineralisation within human tooth cementum identified by secondary ion mass spectrometry[END_REF][START_REF] Mani-Caplazi | Imaging of the human tooth cementum ultrastructure of archeological teeth, using hard x-ray microtomography to determine age-at-death and stress periods[END_REF], 2019;[START_REF] Newham | Exploring the use of tomography for the quantification of cementum growth patterns across the mammal phylogeny[END_REF]. However, if we want to build a robust method, we need to proceed with caution. Cementum studies history (Chapter 1) should be a vivid reminder that rushed publications with oversimplified conclusions can confuse the wider anthropological community not familiar with the intricacies of cementochronology. For example, because aging processes are nonlinear and subject to local influences, using a simple growth ratio based on undefined layers in cementum seems hazardous at best [START_REF] Cerrito | Parturitions, menopause and other physiological stressors are recorded in dental cementum microstructure[END_REF]). Our view is that a proper cementum-specific protocol should be used to test external influences on AEFC growth variation. On the other hand, recent efforts by Mani-Caplazi et al. ( 2019 ) demonstrate a desire to build a robust protocol step by step. This is the way. We hope that these combined efforts will lead to a better understanding of cementum variations.

Noninvasive Virtual Cementochronology

Groundbreaking improvements have been made using synchrotron microCT 3D reconstruction of cementum (Chapters 15 and 16). The possibility of precise and accurate noninvasive analyses, especially when fully automated [START_REF] Newham | Reptile-like physiology in Early  Jurassic stem-mammals[END_REF]Chapter 16), opens a whole new research area for paleoanthropologists and paleontologists who can now extend their questions to the fossil record confidently. Even though synchrotron access remains exceptional, results from such analyses elegantly validate several hypotheses with the most precise and accurate imaging technique available today.

Automated Counting of Increments

Automated counting software is the next area where increased research will provide technological improvements to cementochronology for serial analyses. The development of an automated software has been long in the making (see discussion in Chapter 1), but it is now available for 3D imaging data [START_REF] Newham | Reptile-like physiology in Early  Jurassic stem-mammals[END_REF]Chapter 16). The University of Lille and Leica have been developing a commercial software [START_REF] Bertrand | Computerized cementochronology taking the 16-bit between the teeth[END_REF], which is unfortunately still untested and not publicly available. Conversely, Newham and Naji are currently testing the translation of Newham's 3D algorithm for 2D photomicrographs analysis in a freeware package.

Protocols -The Rise of Cementochronology

With recent methodological improvements and the rise of automation in image analysis, identifying and counting incremental layers in semi-thin/ground sections is now a reliable option for serial analysis of large collections with species-specific standardized protocols. Furthermore, the groundbreaking evolution of virtual cementochronology opens new opportunities to study the fossil record noninvasively. Applications of these new sets of tools on past hominid populations and other fossils open some of the most exciting research opportunities in cementochronology today. However, we must emphasize that these protocols require a significant learning curve to be implemented consistently. Today, few anthropology labs are adequately equipped, and the adaptation of a standardized protocol to the available equipment is counterproductive and should be avoided. By definition, changing any parameter no longer reflects the published protocol but creates a custom version that cannot guarantee tested outcomes.

Part III: Applications

The breadth of applications across disciplines is a testament to the potential of cementochronology as a tool to investigate age and season at death, seasonal mobility patterns, and eventually life-history events.

Seasonality: Discussing Mobility Pattern

In zooarchaeology, cementum increment studies have primarily focused on establishing the season of death of animal dental remains. One of the first major regional studies using cementum that significantly shifted mobility pattern interpretation discussed the seasonal hunting of reindeer during the Upper Paleolithic in France [START_REF] Gordon | Of Men and Reindeer Herds in French Magdalenian Prehistory[END_REF]. Since then, cementochronology has mainly been used to discuss animal migrations (Burke 1993), hunting behavior (Pike-Tay 1991), and mobility patterns [START_REF] Lieberman | Variability in hunter-gatherer seasonal mobility in the Southern Levant: From the Mousterian to the Natufian[END_REF]. With the method's generalization, the studies have progressively expanded to answer more specific questions such as the site function or the specific hunting strategies developed by past human populations (Chapter 17). Others (Chapter 18) propose a regional paleoecological approach to reconstruct past animal coevolution with humans during the Paleolithic. The next stage in cementochronological research applied to zooarchaeologicalcontexts will probably include multiproxies' analyses. Recently Sanchez-Hernandez and colleagues proposed a cross-study involving cementochronology and micro-wear analyses to discuss the site's seasonality of use and the environmental condition during its occupations [START_REF] Sanchez-Hernandez | Combined dental wear and cementum analyses in ungulates reveal the seasonality of Neanderthal occupations in Covalejos Cave (Northern Iberia)[END_REF]. Several authors of this volume obtained a sizable French grant (the ANR DeerPal project, E. Discamps) to characterize past reindeer behavior through a multiproxy analysis including isotopes, cementum, and micro-wear. We expect similar projects to start with the development of new analytical methods that allow us to directly sample elemental signatures of individual cementum increments to characterize the lives of past mammals' populations.

Life-History Events and Fertility

Many researchers have seen the circannual rhythm and continuous growth of acellular cementum throughout life as offering an exciting window of potential information of life-history recorded in the shape, texture, and chemistry of its increments. A host of studies have thus been presented studying the relationship between these factors and various life-history events affecting mammal physiology, most notably pregnancy and parturition (Chapter 8). Ongoing projects from some of this volume's contributors (Newham, Colard, Stock, Price, and Naji) are currently funded to model past human population's fertility rigorously on modern, dedicated cementum databases.

Life Span and Mammalian Evolution

The development of noninvasive 3D methods for imaging cementum increments using synchrotron radiation sources provides volumetric datasets of submicrometer resolution and a consequent precise age estimate for hominid and mammal fossils (Chapters 15,16). Cementum has been imaged successfully in samples from modern humans to Neanderthals [START_REF] Newham | Exploring the use of tomography for the quantification of cementum growth patterns across the mammal phylogeny[END_REF]. Furthermore, Newham and colleagues demonstrated that AEFC increments are still visible in the 200-millionyear-old teeth of early stemmammals [START_REF] Newham | Reptile-like physiology in Early  Jurassic stem-mammals[END_REF]. This unprecedented discovery about cementum preservation opens opportunities to address evolutionary questions in paleontology, such as the evolution of dental attachment and replacement among early mammals. Similarly, for our Homo lineage, the senescence rate and specifically the rate of osteological aging and hominid life spans can be tested directly. For instance, despite increasing knowledge of the anatomy, ecology, and evolution of early mammals, contradictory evidence limits our understanding of their physiology and the timing in which mammalian endothermy arose. Due to the relationship between longevity and basal metabolic rate (BMR) in living mammals, age estimation using virtual cementochronology can be used as a proxy for BMR in fossil taxa. When included in a multiproxy analysis on two early Jurassic stem-mammals, cementochronology has suggested that they had maximum life spans considerably longer than comparably sized living terrestrial mammals and more similar to those of reptiles. The use of this data to estimate BMR suggests that early stem-mammals had yet to develop the elevated endothermic metabolisms of living mammals [START_REF] Newham | Reptile-like physiology in Early  Jurassic stem-mammals[END_REF]. Additionally, new research perspectives in paleoecology are now foreseeable, such as reconstructing animal populations' feeding habits and mobility by analyzing trace elements and stable isotopes at the seasonal cementum increment scale. Finally, validation of the relationship between increment morphology and discrete lifehistory variables and events in living mammals has the potential to foster investigation of these variables in fossil taxa, promoting a new understanding of how they evolved through deep time.

The Elusive Search for the Ideal Skeletal Age Estimator One of the research program's starting points leading to this collaborative volume was the seemingly simple question of estimating the age at death of a human adult skeleton reliably. A slightly more complex follow-up question was how to estimate the age distribution of a skeletal sample accurately. These questions are overlapping, but only in part. A unique answer with a unique age indicator would be possible if such an indicator could provide a continuous, sequential, and unidirectional aging pattern throughout the life span. Preferably, changes should be observed within a short time interval as compared with the total life span. Also, age-related alterations should not be secondary to pathological, metabolic, or nutritional changes. Finally, the indicator should be valid across species -in other words, a precise and accurate estimate with a negligible error rate for all decades of life. In this scenario, the individual age at death could be aggregated into an unbiased age distribution, and our original questions would overlap entirely. Bocquet-Appel and Masset introduced this concept by defining a theoretical minimum requirement for estimating, without considerable risk of error, the individual age of a skeleton and the death structure to which it belongs [START_REF] Bocquet-Appel | Farewell to paleodemography[END_REF]. The correlation between the biological indicator and the age at death must be higher than r = 0.9 (Bocquet-Appel & Masset 1982, 325, table 1). This theoretical minimum is recommended because "the effect of the anthropological reference population on the estimated structure asymptotically decreases as the correlation between age and biological indicator converges towards one" (Bocquet-Appel & Masset 1982, 326). This strong correlation between a biological feature and age is a critical point. Unfortunately, most commonly employed anthropology methods are based on an osteological/dental age indicator staging system where the stages serve as proxies for age [START_REF] Hoppa | The Rostock Manifesto for paleodemography: The way from stage to age[END_REF]. In other words, the determination of biological age is inferred from variables that are poorly correlated with chronological aging. The reason is that skeletal indicators are based on degenerative, wear-and-tear, senescence mechanisms. These changes are primarily correlated with biomechanical loading (lifestyle, fitness) and environmental factors (nutrition, health). However, they are only loosely correlated to an intrinsic aging mechanism (genetic) or a continuous growth process because developmental ontogeny stops in the third decade of life. This methodological distinction is central to understanding some of the underlying mechanisms distinguishing chronological and biological age. If we step out of our discipline and consider the broader field of aging in biology, various biological clocks have been developed over the past decade to predict human health and longevity by providing a quantitative measure of the rate of aging of an individual and their overall resilience more accurately than chronological age. These clocks include DNA methylation patterns, gene expression patterns, serum proteins, and insulin growth factors glycans [START_REF] Jylhava | Biological age predictors[END_REF]. Recent developments in gerontology are starting to identify exact age estimators from living tissue samples. Some are already percolating in forensic anthropology [START_REF] Lee | Forensic DNA methylation profiling from evidence material for investigative leads[END_REF]. For example, the "Horvath clock" is becoming a benchmark in biomedical studies focused on aging mechanisms [START_REF] Horvath | DNA methylation-based biomarkers and the epigenetic clock theory of ageing[END_REF]. More recently, plasma protein age estimation [START_REF] Lehallier | Undulating changes in human plasma proteome profiles across the lifespan[END_REF]) also yielded precision levels above r2 = 0.95, for all decades of life, with a small error rate. Two of these clock calculators [START_REF] Hannum | Genome-wide methylation profiles reveal quantitative views of human aging rates[END_REF][START_REF] Horvath | DNA methylation age of human tissues and cell types[END_REF]) are currently among the most robust predictors of chronological age. Both show high age correlations (r = 0.96 for Horvath and r = 0.91 for Hannum) and small mean deviations from calendar age (3.6 and 4.9 years, respectively) calculated from large validation samples (n = 8000 for Horvath and n = 656 for Hannum). They also cover the entire adult life span and diverse populations [START_REF] Jylhava | Biological age predictors[END_REF]. These are very promising avenues of research on living tissues. Unfortunately, in forensic anthropology and bioarchaeology, we can only rely on techniques for individual skeletal age applicable to dry osteological or dental remains, with two options: (1) Use a single, standard skeletal/dental age indicator [START_REF] Buikstra | Standards for Data Collection from Human Skeletal Remains: Proceedings of a Seminar at the Field Museum of Natural History, Organized[END_REF]. The resulting estimation is invariably associated with a significant standard error (>10 years at best), especially for adults above fifty. (2) Use a multifactorial method combining several skeletal indicators into an "expert opinion" [START_REF] Milner | Transition analysis: A validation study with known-age modern American skeletons[END_REF]. Unfortunately, few procedures provide a viable mechanism for combining age-related information from separate anatomical structures in different parts of the skeleton in terms of point estimate or confidence intervals. From a statistical standpoint, this is problematic because the information from multiple age indicators is not independent. Without an unrealistically large reference sample, a reliable estimate is impossible, at least not without considerable simplifications [START_REF] Boldsen | Transition analysis: A new method for estimating age from skeletons[END_REF]. However, the transition analysis method's recent evolution (http://ta3info.com) explores an entirely new approach to the expert-opinion aging process. We are looking forward to Milner and Boldsen group's fully published method and the subsequent validation tests. Until then, these two options are highly unsatisfactory, especially with the growing rigorous legal criteria that expert witnesses must uphold in forensic contexts. Recommendations are starting to emerge favoring a probabilistic summary of results, presenting posterior probabilities for each age [START_REF] Chariot | Age estimation in undocumented migrant adolescents: Medical response to judicial authorities[END_REF][START_REF] Konigsberg | Optimal trait scoring for age estimation[END_REF]. We argue that this is the most rigorous approach when using standard anthropological age indicators. In paleodemography, the debate to estimate a skeletal sample age-at-death distribution has come a long way since its initial deconstruction (Bocquet-Appel & Massat 1982). A strong consensus has emerged promoting a Bayesian approach of collective age estimation [START_REF] Bocquet-Appel | Estimation of an age distribution with its confidence intervals using an iterative bayesian procedure and a bootstrap sampling approach[END_REF][START_REF] Hoppa | The Rostock Manifesto for paleodemography: The way from stage to age[END_REF][START_REF] Konigsberg | Multivariate cumulative probit for age estimation using ordinal categorical data[END_REF][START_REF] Seguy | Estimating the age structure of a buried adult population: A new statistical approach applied to archaeological digs in France[END_REF]. In this probabilistic framework, two main trends have emerged. The English-speaking trend is mostly focused on parametric hazard models [START_REF] Gage | Variation and classification of human age patterns of mortality: Analysis using competing hazards models[END_REF][START_REF] Konigsberg | Optimal trait scoring for age estimation[END_REF][START_REF] Usher | Reference samples: The first step in linking biology and age in the human skeleton[END_REF]Chapter 22). In contrast, the French school is more data driven and proposes a fully Bayesian algorithm [START_REF] Caussinus | Un probleme cle de la paleodemographie: Comment estimer l'age au deces?[END_REF]Caussinus & Courgeau 2013). This method integrates any standard skeletal indicator into an age distribution with relatively tight credible intervals. Unfortunately, even if this should not be an excuse, both options' level of statistical proficiency prevents a broader deployment of these methods.

Back to the Root

Of course, this volume offers a third alternative. Cementochronology is, by far, the technique that complies with the most criteria defining a robust age estimation method (Kemkes-Grottenthaler 2002). Cementochronology is (1) based on a continuous growth process; (2) the most independently validated indicator, not only in humans but across eighty different species; (3) the most precise and accurate; (4) independent from any reference collection; and (5) based on teeth, one of the most available remains in archaeological samples. Cementochronology, with its high correlation to chronologicalage (>0.9, see Chapter [START_REF] Bertrand | Age at death estimation by cementochronology: Too precise to be true or too precise to be accurate?[END_REF]), cementochronology has never been described as a "perfect" age indicator. However, several authors have argued (Colard et al. 2015;Naji et al. 2016b;Introduction) that there is a clear distinction between cementochronology, which measures an actual growth process, and other skeletal adult age indicators, all based on degenerative changes. Cementum growth is subject to some aging processes and individual variations of phenotypic expression. Cementochronology estimates should thus always be assorted with an error rate. However, to date, this error rate is by far the smallest of any other skeletal or dental age indicators for any adult decade of life across species (Chapter 1: Table 1.1). Surely, this point is worth emphasizing. Theoretically, it is probably unreasonable to believe that a perfect indicator exists considering the natural variation expressed in any biological organism's growth/ senescent processes. In particular, recent discoveries on aging mechanisms [START_REF] Farr | The spectrum of fundamental basic science discoveries contributing to organismal aging[END_REF][START_REF] Lopez-Otin | The hallmarks of aging[END_REF] suggest that none of the visible changes associated with age is a linear process. Maintenance of body function and cohesion is contextdependent and subject-dependent, predicting a nonlinear relationship between any biological and chronological age, which we have all empirically observed [START_REF] Jylhava | Biological age predictors[END_REF][START_REF] Woodward | Skeletochronology[END_REF]). Any such pattern can and should be integrated into age estimation modeling to increase the method's precision, especially in the last decades of our life span.

Cementochronology, the Final Frontier

Traditionally, one theoretical and two practical concerns were usually invoked to dismiss cementochronology:

(1) the evolutionary process explaining circannual growth cycles, (2) the subjectivity of counting the increments, and (3) the destructive nature of histological thin-section preparation. We argue that the broader chronobiology framework (Introduction, Part I) provides the necessary hypotheses to test the first theoretical point further. Moreover, with the recent development of species-specific standardized protocols, the development of counting automation, and 3D noninvasive analysis (Part II), we argue that the two primary practical sources of concerns (subjectivity of counts and destruction of samples) have viable solutions. With these recent advances, we argue that cementochronology has the potential to renew several anthropology themes robustly (Part III). The goal of this volume was to promote interdisciplinary discussions and stimulate new research opportunities in cementochronology. The enthusiastic collaborations resulting from this endeavor illustrate the success of the project. We hope that the Cementochronology Research Program, and anyone interested in cementum, will continue their collaborative effort and interdisciplinary approach to expand the boundaries of cementum in anthropology.

  ), Colard et al. (2015) (Chapter 14), and Newham et al. (2020) (Chapters 8 and 16

  1, Table 1.1), should logically expand as a viable age indicator to estimate individual age at death or collective age distribution in all contexts (Chapters 19-23). Nevertheless, the search for an ideal skeletal age indicator (progressive, precise, accurate, universal, and simple) continues, motivated by the multitude of applications in wildlife biology, forensic anthropology, bioarchaeology, paleodemography, paleopathology, zooarchaeology, and paleontology. Contrary to what some authors claim

Acknowledgments

We would like to thank the ANR CemeNTAA (ANR-14-CE31-0011) and the ANR DeerPal (ANR-18-CE03-0007) for funding some of the experiments, travels, symposiums, and articles that led to the publication of this edited volume. We thank Elis Newham for help on an earlier draft. We would like to thank all the contributors who suffered our relentless reminders. We would also like to thank Daniel Antoine for his personal involvement in motivating us to seek a publisher during the 2017 AAPA symposium.


Marc Baaden

Mark S P Sansom
email: mark.sansom@biop.ox.ac.uk

OmpT: Molecular Dynamics Simulations of an Outer Membrane Enzyme

Five MD simulations (total duration >25 ns) have been performed on the E. coli outer membrane protease OmpT embedded in a DMPC lipid bilayer. Globally the protein is conformationally stable. Some degree of tilt of the b-barrel is observed relative to the bilayer plane. The greatest degree of conformational flexibility is seen in the extracellular loops. A complex network of fluctuating H-bonds is formed between the active site residues, such that the Asp210-His212 interaction is maintained throughout, whilst His212 and Asp83 are often bridged by a water molecule. This supports a catalytic mechanism whereby Asp83 and His212 bind a water molecule which attacks the peptide carbonyl. A configuration yielded by docking calculations of OmpT simulation snapshots and a model substrate peptide Ala-Arg-Arg-Ala was used as the starting point for an extended Hückel calculation on the docked peptide. These placed the LUMO mainly on the carbon atom of the central C=O in the scissile peptide bond, thus favouring attack on the central peptide by the water held by residues Asp83 and His212. The trajectories of water molecules reveal exchange of waters between the intracellular face of the membrane and the interior of the barrel but no exchange at the extracellular mouth. This suggests that the 'pore-like' region in the centre of OmpT may enable access of water to the active site 'from below'. The simulations appear to reveal the presence of specific lipid interaction sites on the surface of the OmpT barrel. This reveals the ability of extended MD simulations to provide meaningful information on protein/lipid interactions.

Introduction

Bacterial outer membrane proteins (OMPs) all share a common architecture, that of a transmembrane domain formed by an anti-parallel b-barrel. Up to date the structures of 20 such OMPs have been solved by X-ray diffraction [START_REF] Buchanan | ß-Barrel proteins from bacterial outer membranes: structure, function and refolding[END_REF][START_REF] Koebnik | Structure and function of bacterial outer membrane proteins: barrels in a nutshell[END_REF] and by NMR [START_REF] Arora | Structure of outer membrane protein A transmembrane domain by NMR spectroscopy[END_REF][START_REF] Fernandez | Solution NMR studies of the integral membrane proteins OmpX and OmpA from Escherichia coli[END_REF]. Thus, OMPs provide an opportunity for molecular dynamics (MD) simulation studies to explore the conformational dynamics of a whole family of structurally related membrane proteins, in order to define both common dynamic properties and functionally important differences between individual species of OMPs (Domene et al., 2003a;[START_REF] Bond | The simulation approach to bacterial outer membrane proteins[END_REF].

Perhaps the best characterised family of OMPs are the porins [START_REF] Cowan | Bacterial porins -lessons from 3 high-resolution structures[END_REF][START_REF] Schirmer | General and specific porins from bacterial outer membranes[END_REF][START_REF] Achouak | Multiple facets of bacterial porins[END_REF]. These include both relatively non-specific general diffusion pores and also more specific passive pores (e.g. for oligosaccharides [START_REF] Schirmer | Structural basis for sugar translocation through maltoporin channels at 3.1A resolution[END_REF][START_REF] Forst | Structure of the sucrose-specific porin ScrY from Salmonella typhimurium and its complex with sucrose[END_REF]) across the outer membrane. Thus, they are an important component of the transport properties of the bacterial membrane. Other transport proteins in outer membranes include those for ferric ions [START_REF] Locher | Transmembrane signalling across the ligand-gated FhuA receptor; crystal structures of free and ferrichrome-bound states reveal allosteric changes[END_REF][START_REF] Ferguson | Siderophoremediated iron transport: crystal structure of FhuA with bound lipopolysaccharide[END_REF][START_REF] Buchanan | Crystal structure of the outer membrane active transporter FepA from Escherichia coli[END_REF] and for vitamin B12 [START_REF] Chimento | Substrate-induced transmembrane signaling in the cobalamin transporter BtuB[END_REF], and export pathways for polypeptide toxins and hydrophobic drugs such as TolC [START_REF] Koronakis | Crystal structure of the bacterial membrane protein TolC central to multidrug efflux and protein export[END_REF].

In addition to transport proteins, outer membranes include a number of membrane bound enzymes. Several structures of such enzymes have been determined, including those of a protease OmpT [START_REF] Vandeputte-Rutten | Crystal structure of the outer membrane protease OmpT from Escherichia coli suggests a novel catalytic site[END_REF], and two OMPs acting on lipid substrates, OMPLA [START_REF] Snijder | Structural evidence for dimerization-regulated activation of an integral membrane phospholipase[END_REF] and PagP [START_REF] Hwang | Solution structure and dynamics of the outer membrane enzyme PagP by NMR[END_REF].

In the present work we focus on OmpT, the Escherichia coli outer membrane endoprotease Omptin (EC 3.4.21.87), which shows maximum enzymatic activity for cleavage sites with two consecutive basic amino acids (Arg-Arg, Lys-Arg, Lys-Lys) [START_REF] Sugimura | Purification, characterization, and primary structure of Escherichia coli protease VII with specificity for paired basic residues: identity of protease VII and OmpT[END_REF]. OmpT also cleaves a number of more remote sequences as has been shown by activity assays and library screening [START_REF] Dekker | Substrate specificity of the integral membrane protease OmpT determined by spatially addressed peptide libraries[END_REF]. It is functional as a monomer and has an autoproteolytic site at K217-R218. Mutation of the three residues S99, G216 and K217 was a prerequisite for the crystallographic structure determination [START_REF] Vandeputte-Rutten | Crystal structure of the outer membrane protease OmpT from Escherichia coli suggests a novel catalytic site[END_REF].

The OmpT crystal structure comprises 297 amino acids and has revealed its prolonged 10stranded beta-barrel architecture with a central elliptical crosssection of ca. 13 x 16 Å. The active site sits in the extracellular half of the structure at the base of two long loops that protrude from the barrel flanking a central binding pocket. Four residues have been identified as essential part of the active site, D210/H212 on one side of the pocket and D83/D85 on the opposite side.

The biological function of OmpT is not fully established, but it may have a protective role in pathogenic E. coli. OmpT degrades a variety of positively charged antimicrobial peptides and might be involved in urinary tract disease and DNA excision repair.

MD simulations provide a valuable tool for studying membrane proteins, enabling us to probe their conformational dynamics in both membrane and detergent micelle environments [START_REF] Bond | Membrane protein dynamics vs. environment: simulations of OmpA in a micelle and in a bilayer[END_REF]. They are of particular value in enabling us to extrapolate from the essentially static (time and space averaged) structure revealed by X-ray diffraction to a more dynamic picture of the behaviour of a single OMP molecule in a more realistic environment mimicking a small patch of the bacterial membrane. MD simulations have been employed in a number of studies of OMPs, most notably to probe protein and solvent dynamics in relationship to permeation mechanisms in porins [START_REF] Tieleman | A molecular dynamics study of the pores formed by Escherichia coli OmpF porin in a fully hydrated palmitoyloleoylphosphatidylcholine bilayer[END_REF][START_REF] Im | Ions and counterions in a biological channel: A molecular dynamics simulation of OmpF porin from Escherichia coli in an explicit membrane with 1 M KCl aqueous salt solution[END_REF], to explore possible pore gating mechanisms in OmpA [START_REF] Bond | OmpA -A pore or not a pore? Simulation and modelling studies[END_REF][START_REF] Bond | Membrane protein dynamics vs. environment: simulations of OmpA in a micelle and in a bilayer[END_REF], to explore dynamics in relationship to transport in FhuA [START_REF] Faraldo-Gómez | Molecular dynamics simulations of the bacterial outer membrane protein FhuA: a comparative study of the ferrichrome-free and bound states[END_REF], and to examine the role of calcium binding and dimerisation in the catalytic mechanism of OMPLA [START_REF] Baaden | A molecular dynamics investigation of mono-and dimeric states of the outer membrane enzyme OMPLA[END_REF].

In the current study we employ MD simulations to examine the conformational dynamics of the active site of the outer membrane protease OmpT. We also explore the interactions of OmpT with phospholipid molecules, which is of some interest in the context of possible specific interactions of OmpT with lipid A and their role in stability and/or function of OmpT.

Methods

System Preparation

Repair of protein structures

We started from chain A of the 2.6 Å crystal structure (PDB code 1i78) where missing sidechain atoms were modelled using a rotamer library [START_REF] Chinea | The use of position specific rotamers in model building by homology[END_REF]. We reversed the mutations required for the crystallization of OmpT: A99S, K216G and G217K.

pKA Calculations

It has been our experience 66 that more stable MD simulations can be obtained if pKA calculations are used to aid assignment of ionisation states to acidic and basic residues in a protein. Thus, pKA calculations for OmpT were performed using the program WhatIf 67 . This approach combines calculation of the energies of different protonation states of ionisable residues via the Poisson-Boltzmann equation with modelling of local structural changes associated with (de)protonation of ionisable residues [START_REF] Nielsen | Optimizing the hydrogen-bond network in Poisson-Boltzmann equation-based pK(a) calculations[END_REF]. To mimic the dielectric environment presented by a membrane pKA calculations were performed on OmpT embedded in a slab of methane molecules. The result was a series of titration curves for each ionisable residue, from which protonation states were assigned.

Insertion into a Bilayer

OmpT was inserted into a pre-equilibrated, hydrated (65 waters per lipid molecule) DMPC bilayer using the methods described by [START_REF] Faraldo-Gómez | Setup and optimisation of membrane protein simulations[END_REF]. Briefly, a hole of approximately the same size as the protein is made in the lipid bilayer by removal of DMPC molecules. A simulation is then run in which the lipid molecules at the edge of the hole experience local forces perpendicular to a molecular surface generated from the protein to be inserted. This allows the lipids to relax around a protein-shaped cavity into which the protein is then introduced. Na + and Cl -ions equivalent to a concentration of ca. 0.1 M were then added. A typical view of the total simulation system is shown in Fig. 1C.

Simulations

All simulations were performed using GROMACS 3.1.4. [START_REF] Berendsen | GROMACS: A messagepassing parallel molecular dynamics implementation[END_REF] The force field used was GROMOS-87. (van Gunsteren and Berendsen, 1987) The SPC water model was employed. [START_REF] Berendsen | Intermolecular Forces[END_REF] The simulations were run in the NPT ensemble. The temperature was maintained at 310 K using a Berendsen thermostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] with τ = 0.1 ps. The pressure was maintained at 1 bar by anisotropically coupling x, y and z components to a Berendsen barostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] with τ = 1.0 ps and compressibility of 4.5 10 -5 bar --1 in all three dimensions. The time step for integration was 2 fs, and coordinates were saved every 1 ps for subsequent analysis. Electrostatic interactions were calculated using either: a cutoff of 18 Å and van der Waals interactions were truncated at 14 Å; or with PME [START_REF] Darden | Particle mesh Ewald -an N.log(N) method for Ewald sums in large systems[END_REF]) using a 10 Å cutoff for both the real-space calculation and the van der Waals interactions. The LINCS algorithm [START_REF] Hess | LINCS: A linear constraint solver for molecular simulations[END_REF] was used to constrain all bond lengths.

Following insertion of the protein into the bilayer (see above) an equilibration simulation was performed during which restraints on the protein atoms were gradually released. This equilibration period was 0.6 ns. Production simulations of up to 10 ns duration were then performed.

Peptide Docking

Docking calculations of the three peptides ARRA, AKKA and AK(D)RA were performed with autodock3 [START_REF] Morris | Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function[END_REF] for seven protein conformations, identified by cluster analysis of all trajectories. In a second round we focused on one of the protein conformations, where protonation states were changed to better reflect a potential transition state, Asp83 being protonated and His212 neutral. A water molecule was explicitly added between Asp83 and His212. We selected one of the ARRA complexes and refined sidechain positions interactively with the Yasara software [START_REF] Krieger | Increasing the precision of comparative models with YASARA NOVA: a self-parameterizing force field[END_REF]. Protein-substrate interactions were analysed with the program LIGPLOT [START_REF] Wallace | LIGPLOT: a program to generate schematic diagrams of protein-ligand interactions[END_REF]. The STC software [START_REF] Lavigne | Structure-based thermodynamic analysis of the dissociation of protein phosphatase-1 catalytic subunit and microcystin-LR docked complexes[END_REF] was used to estimate the free energies of active site/peptide interactions via a structurebased thermodynamic approach. Given the simple docking procedure used, the approximate nature of the resulting complex and the short length of the peptide, we only derived relative contributions of protein sidechains to the binding energy.

Additional Analysis Programs

Secondary structure elements were identified using the program DSSP [START_REF] Kabsch | Dictionary of protein secondary structure:-pattern-recognition of hydrogen-bonded and geometrical features[END_REF], graphical representations were prepared with VMD [START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF] and Raster3D [START_REF] Merritt | Raster3D: Photorealistic molecular graphics[END_REF]. Water bridges were identified with an adapted Gromacs analysis program provided by JD Faraldo-Gómez. Lipid-protein interactions were analysed using the Proximus database and associated tools (Deol and Sansom, unpublished results).

Results

Simulations and Stability

The simulation system, consisting of an OmpT monomer embedded in a DMPC bilayer, is shown in Fig. 1. In determining the initial configuration of this system, attention was paid to the location of the bands of amphipathic aromatic amino acid (i.e. Trp and Tyr) sidechains on the surface of the b-barrel, as it is thought that such residues help to 'lock' integral membrane proteins within a lipid bilayer [START_REF] Killian | How proteins adapt to a membrane-water interface[END_REF][START_REF] Lee | Lipid-protein interactions in biological membranes: a structural perspective[END_REF][START_REF] Fyfe | Probing the interface between membrane proteins and membrane lipids by X-ray crystallography[END_REF]. The other class of residues that have been suggested to interact with lipid headgroups are basic sidechains. In this context it is noteworthy that there are many basic residues in the upper (i.e. extracellular) half of the OmpT barrel, pointing out towards the lipid headgroups (Fig. 1A: also see discussion below).

Five simulations were performed (Table 1). OMPT1 and OMPT2 were both of 10 ns duration, but differed in the manner in which long-range electrostatic interactions were approximated. OMPT3a-c were short (2.5 ns) reruns of the OMPT2 simulation using different random seeds and starting conformations. By performing multiple simulations we aimed to explore the robustness of our results to changes in simulation protocol (OMPT1 vs. OMPT2) and to stochastic fluctuations and differences between multiple runs of a given simulation protocol.

Furthermore, better sampling is achieved by multiple simulation runs [START_REF] Caves | Locally accessible conformations of proteins: Multiple molecular dynamics simulations of crambin[END_REF].

Visualisation of the results of e.g. simulation OMPT2 (Fig. 1C) suggests that whilst there is no overall translation of the OmpT b-barrel relative to the lipid bilayer, some degree of tilting occurs. This is discussed in more detail below. However, extended simulations of phospholipid bilayers have revealed surface undulations [START_REF] Lindahl | Mesoscopic undulations and thickness fluctuations in lipid bilayers from molecular dynamics simulations[END_REF] and so we cannot exclude the possibility of similar long (i.e. >> 10 ns) timescale fluctuations in OMP orientations in membranes.

The overall drift of OmpT from the initial (i.e. crystallographic) structure provides a measure of the conformational stability of this protein in a membrane environment. Drift can be simply measured in terms of the Cα root mean square deviation (RMSD) from the starting structure. In both simulations OMPT1 and OMPT2 (Fig. 2) it can be seen that the b-barrel is very stable, with a final Cα RMSD of ~1 Å, whereas the highest degree of structural drift is seen in the loop regions, for which the cutoff simulation (i.e. OMPT2) may show a slightly higher degree of drift on a 10 ns timescale.

Detailed analysis of the convergence of the simulations using principle components analysis (PCA) and other methods will be published elsewhere [START_REF] Faraldo-Gómez | Conformational sampling and dynamics of membrane proteins from 10-nanosecond computer simulations[END_REF]. To summarize, the OmpT simulations were among the best converged systems comparing 9 different membrane proteins in a lipid environment. The similarity of both halves of the simulation -a measure for the sampling completeness -derived from PCA is 81 % for the transmembrane part and 66 % for the whole protein. This compares very favourably to shorter simulations of the OmpF porin with values from 12 to 16 % [START_REF] Hess | Convergence of sampling in protein simulations[END_REF]. The structural drift of the backbone of the four active site residues was found to be less than 1 Å.

Analysis of the radius of gyration confirms that OmpT is a stable protein in multinanosecond simulations and suggests some degree of tilting of OmpT in the bilayer, which is more pronounced with OMPT1. More exact protein tilt analysis confirms the significant tilting of OmpT in the bilayer (this is discussed below). Overall, the simulations confirm that, despite the low resolution (2.6 Å), the OmpT structure is remarkably stable in simulations on a ~10 ns timescale.

Conformational Fluctuations

Whilst the OmpT structure is globally stable in the simulations, this is not to imply that no significant conformational fluctuations take place. Thus, it is of interest to also examine the magnitude of the conformational fluctuations in different regions of the structure and to compare these with the experimental B-values (whilst not forgetting the resolution of the structure, namely 2.6 Å). In Fig. 3A the Cα root mean square fluctuations (RMSFs) are shown as a function of residue number for OMPT1, and compared with the equivalent RMSFs derived from the B-values. Qualitatively the two curves agree, with the highest fluctuations being seen in the extracellular loops and the lowest fluctuations in the core of the b-strands. Two major differences are observed. Firstly the peak values of the RMSFs for the loops are higher in the simulation than in the crystal structure. This may reflect constraints on loop mobility present in the crystal which are removed when the protein is in a bilayer. Such an interpretation is supported by recent simulations of the small outer membrane protein OmpA in a crystal vs. a bilayer environment (Bond and Sansom, ms. in preparation). The second major difference is that the RMSFs derived from the B-values are significantly higher than those seen in the simulations for the core regions of the b-barrel. This probably reflects some contribution of static crystal disorder (e.g. mosaicity)

to the B-values as well as undersampling of long-timescale motions of the barrel [START_REF] Faraldo-Gómez | Conformational sampling and dynamics of membrane proteins from 10-nanosecond computer simulations[END_REF].

We have also examined the Cα RMSF values as a function of the z coordinate of the Cα atoms (i.e. the position of the Cα atom projected onto the approximate bilayer normal). For all three simulations this analysis (data not shown) reveals a similar pattern to that seen in other outer membrane protein simulations [START_REF] Bond | The simulation approach to bacterial outer membrane proteins[END_REF], namely the RMSFs are lowest (~0.6 Å) in the centre of the membrane (also the middle of the barrel) and rise at either end, being highest in the extramembraneous loops.

Visualisation (Fig. 3B) of superimposed snapshots taken throughout the simulations suggests that the pattern of mobility defines two relatively rigid 'halves' of the barrel, with highly mobile extracellular loops (two of which contain the active site residues sitting at the region connecting with the barrel) flanking a potential peptide binding site in between. This poses the question of the possible role(s) of such loop flexibility in the catalytic mechanism of OmpT. One possibility is that the mobile loops enhance (i.e. lower the activation energy for) entry/exit of substrate/product from the active site.

We have also monitored the secondary structure as a function of time for each simulation. Comparing the three simulations (i.e. treating the OMPT3a-c simulations together) and examining the predominant secondary structure for each residue during each simulation (Fig. 3C), one can see that differences in secondary structure between the three simulations are only observed for the extracellular loops, in agreement with the pattern of fluctuations described above.

The Active Site

The mechanism proposed for OmpT has active site residues located on two extracellular loops on opposite faces of the barrel (cf. Fig. 1A), namely a His212-Asp210 dyad on one side of the proposed peptide binding site, and an Asp83-Asp85 couple on the other side. The His212-Asp210 dyad is thought to activate a water molecule for nucleophilic attack on the C of the peptide bond whilst the Asp83-Asp85 couple may have a dual role, both by coordinating and thus helping to orient the nucleophilic water [START_REF] Vandeputte-Rutten | Crystal structure of the outer membrane protease OmpT from Escherichia coli suggests a novel catalytic site[END_REF] and possibly by stabilizing the oxyanion intermediate via a shared proton [START_REF] Kramer | Identification of essential acidic residues of outer membrane protease OmpT supports a novel active site[END_REF]. Earlier mutagenesis studies [START_REF] Kramer | Identification of active site serine and histidine residues in Escherichia coli outer membrane protease OmpT[END_REF] had suggested a role for Ser99 although this is less clear in terms of the structure. It has been suggested [START_REF] Kramer | Identification of essential acidic residues of outer membrane protease OmpT supports a novel active site[END_REF]) that Glu27 and Asp208 may also play a role in substrate binding via interaction with a substrate arginine sidechain. The simulations of OmpT provide an opportunity to examine the structural integrity of the active site and the nature of the interactions between the key residues as they evolve over time, and in particular to explore the interactions of water molecules with active site sidechains.

There is a complex network of H-bonds (Fig. 4A/B) between the active site residues, which fluctuates as a function of time (Fig. 4C; Table 2). Firstly, in all of the simulations the Asp210-His212 H-bonding interaction is maintained throughout. His212 and Asp83 are often bridged by a water molecule (Fig. 4) thus supporting the mechanism of [START_REF] Kramer | Identification of essential acidic residues of outer membrane protease OmpT supports a novel active site[END_REF] where Asp83 and His212 are proposed to bind the water which goes on to attack the peptide carbonyl. Occasionally, Asp83 and His212 are bridged by a chain of two water molecules.

Interestingly, the Asp83-water-His212 interaction is sometimes replaced by an Asp83-Ser99 interaction. The loss of enzyme activity on mutating Ser99 to Ala suggests that this alternative Asp-Ser interaction may play a functionally important role. The Asp83-Asp85 pair interact with one another via one or more bridging water molecules. Finally, Glu27 switches between interacting with Asp208 directly and with Asp210 via an intervening water. Thus, water molecules seem to play an important role in the active site, both in terms of being a possible component of the catalytic mechanism per se (the water(s) held by His212 and Asp83) and by bridging between sidechains to maintain the integrity of the active site. Of course, these simulations are in the absence of bound peptide substrate, and a different configuration of sidechains and waters may be expected in the latter state.

As part of the preparation of the simulation system (see Methods) we performed pKA calculations. These suggested that His212 protonation was favoured by ~10% at neutral pH (the calculated pKA was ~7.9). In combination with the results of the H-bonding analysis discussed above, this might suggest that a proton transfer shuttle between water, His212 and Asp83 may play a key role in the catalytic mechanism. Interestingly, the pKA calculations also suggested that the pKA of Asp85 was perturbed from its standard value (to ~6.6).

The simulations were all performed in the absence of bound peptide substrate. In order to explore the influence of active site fluctuations on enzyme-substrate interactions we have performed docking simulations (see Methods section) between selected OmpT snapshots from the simulations and three model peptides: Ala-Lys-Lys-Ala, Ala-Arg-Arg-Ala, and Ala-Lys-(D)Arg-Ala. The first two are models of substrates, whereas the third peptide is an inhibitor. The results from docking these peptides to different snapshots of OmpT from the simulations provide a range of complexes. Some of the observed configurations are compatible with the catalytic mechanisms discussed and place the peptide cleavage site near the catalytic residues. However, the precision of the docking calculation does not permit a detailed insight into the mechanism as several alternative docked structures are close to one another in energy. In Fig. 5 we show an interactively refined docked complex of Ala-Arg-Arg-Ala with the active site configuration from simulation OMPT1 at 0 ps. It can be seen that the scissile peptide bond is located between the Asp83-Asp85 and Asp210-His212 pairs. This conformation was chosen, because a water molecule was present between His212 and Asp83. The water is supposed to perform the nucleophilic attack. In order to take into account the effect of a second water between Asp83 and Asp85, we protonated Asp83, which is therefore able to provide a proton for the stabilization of the oxyanion produced during hydrolysis. The N-terminal Arg of the peptide interacts with Glu27 (distance 2.8 Å) and Asp208 (distance 3.0 Å) and the C-terminal Arg with Asp97 (distance 3.1 Å). Furthermore the N-terminal Ala interacts with Met81 and Ile170 as has been hypothesized previously [START_REF] Vandeputte-Rutten | Crystal structure of the outer membrane protease OmpT from Escherichia coli suggests a novel catalytic site[END_REF]. Thus the results of the docking are consistent with the mechanism proposed on the basis of the crystal structure and mutagenesis studies.

Another interesting feature of the OmpT active site are the fluctuations of the putative peptide binding cleft, in particular the hot spot delineated by the four residues D83, D85, N210 and H212. One flexibility measure is the evolution of the H212-D83 and H212-D85 distances vs. time, which fluctuate between 3 and 8 Å (graphs not shown) and show a clear bimodal distribution with maxima at 4.0 and 5.6 Å in simulation OMPT2.

To further examine the peptide docked at the active site we performed an extended Hückel calculation on the Ala-Arg-Arg-Ala peptide in the docked conformation discussed above. The results of this calculation (data not shown) placed the LUMO mainly on the carbon atom of the central C=O in the scissile peptide bond. A control calculation on an energy minimized conformation of the isolated peptide -with no docking constraints -placed the LUMO in a different position incompatible with the expected enzymatic reaction. Thus the docked conformation seems to favour attack on the central peptide bond and the preferred direction of attack would seem to be from "below", i.e. from the direction of the water held by residues Asp83 and His212.

The results of a thermodynamic analysis (described in the Methods section) of the peptide/OmpT complex shown in Fig. 5 suggests that the predominant interactions between the peptide and the binding site include those of Glu27 and Asp210, but also, more surprisingly, those of Tyr150, Arg168, Lys217 and Tyr221.

The results of our simulations seem to be broadly compatible with suggested mechanism(s) [START_REF] Kramer | Identification of essential acidic residues of outer membrane protease OmpT supports a novel active site[END_REF]. In particular, we have evidence for the stability of the Asp210-His212 dyad, and for the presence of one -occasionally two -bridging water molecules held by His212 and Asp83. This strongly suggests that the water molecule held by His212 and activated by the catalytic dyad is 'poised' to attack the scissile peptide bond, the peptide being oriented such that its LUMO is directed towards the water oxygen. As to the stabilizing role of Asp83 and Asp85, several possibilities can be considered. Firstly, the water held between Asp83 and Asp85 may act as the donor of a shared proton between these two residues, in agreement with the mechanism proposed by [START_REF] Kramer | Identification of essential acidic residues of outer membrane protease OmpT supports a novel active site[END_REF]. Secondly, the attacking water molecule could simultaneously provide a proton to stabilize the oxyanion intermediate and hydrolyse the peptide bond. Thirdly, if two water molecules were bridging, one could play a stabilizing role whilst the other acts as a nucleophile.

Water

The crystal structure of OmpT reveals between 6 and 9 waters within the OmpT barrel (depending on which of the two monomers in the asymmetric unit one examines). Given the behaviour of water molecules within the pore-like β-barrel of OmpA we were interested to examine the dynamic properties of water within the OmpT barrel. In particular, we wished to examine whether the water inside the β-barrel behaved simply as a 'structural' element or whether OmpT could form a potential water permeable pore through the outer membrane.

The time average distribution of water molecules inside the barrel (Fig. 6A) indicates that water can access most of the interior of the barrel but that there is a region of low water density around the extracellular mouth of the barrel. This low water density region corresponds to the narrowest region of the barrel in the crystal structure, where the barrel is 'squashed' such as to have a highly elliptical cross-section. If we examine the trajectory of a single (typical) water molecule (Fig. 6B) we can see ready exchange of the water molecule between the intracellular face of the membrane and the interior of the barrel but no exchange of water at the extracellular mouth. This is supported by following the trajectories of a number of water molecules within the barrel projected onto the z (i.e. approximate barrel) axis (Fig. 6C).

What does this mean in terms of the function of OmpT? One is tempted to speculate that the 'pore-like' region in the centre of the OmpT molecule may allow access of water to/from the active site 'from below' i.e. on the side of the active site facing towards the interior of the βbarrel just beneath the putative peptide binding site.

Lipid/Protein Interactions

The OMPT simulations provide an opportunity to explore further the nature of lipidprotein interactions in a lipid bilayer mimicking the bacterial outer membrane and the extent to which these may be characterised by simulations. This is a topic of some general importance, as an improved understanding of protein/lipid interactions may enable us to better predict the structure of membrane proteins. Furthermore, interactions between lipid molecules and membrane proteins are known to play important roles in the stability and structural integrity of membrane proteins [START_REF] Killian | How proteins adapt to a membrane-water interface[END_REF][START_REF] Lee | Lipid-protein interactions in biological membranes: a structural perspective[END_REF]. Furthermore, OmpT requires LPS for its function and on the basis of structural comparisons with FhuA [START_REF] Ferguson | Siderophoremediated iron transport: crystal structure of FhuA with bound lipopolysaccharide[END_REF] it has been suggested that there may be a specific LPS binding site on the surface at the extracellular end of the OmpT β-barrel [START_REF] Vandeputte-Rutten | Crystal structure of the outer membrane protease OmpT from Escherichia coli suggests a novel catalytic site[END_REF].

As discussed briefly above we have observed a degree of tilting of the OmpT barrel relative to the lipid bilayer. For example, at the end of simulation OMPT1 the barrel axis is tilted ~20º relative to the (overall) bilayer plane. We have seen a similar tilting of the barrel of the simple outer membrane protein OmpA [START_REF] Bond | OmpA -A pore or not a pore? Simulation and modelling studies[END_REF] and a similar tilting of OmpA has been predicted by [START_REF] Basyn | Prediction of membrane protein orientation in lipid bilayers: a theoretical approach[END_REF] using a simplified potential for protein/membrane interactions. We therefore suggest that the tilt of OmpT may reflect an underlying asymmetry in the interaction of the outer barrel envelope with the bilayer. This is suggested by the agreement between the tilt seen in MD simulations and that generated using a simple hydrophobicity potential for sidechain/bilayer interactions (Deol & Sansom, unpublished results). It is possible that the development of the tilt of OmpT relative to the bilayer during the simulations may be related to the increase in the number of lipid/protein H-bonds observed during the early stages of both simulations OMPT1 and OMPT2 (Fig. 7).

This increase in number of H-bonds is suggested to correspond to a slow (> 2 ns) optimisation of the interaction of the protein with the bilayer. Similar timescales for fluctuations of number of protein/lipid interactions have been observed in simulations of OmpA and KcsA (Domene et al., 2003b).

As mentioned above, on the basis of comparing the crystal structures of FhuA (with bound LPS; PDB code 1QFG) and of OmpT, [START_REF] Vandeputte-Rutten | Crystal structure of the outer membrane protease OmpT from Escherichia coli suggests a novel catalytic site[END_REF] suggested a LPS binding site for OmpT. This consisted of OmpT residues Tyr134, Glu136, Arg138, Arg175 and Lys226. In the FhuA structure the primary contacts from the protein to the Lipid A portion of LPS are via a lysine and arginine cluster to the diphosphate moiety and via a lysine and glutamine cluster to the single phosphate (Fig. 8A).

We have analysed the principal contacts between the phosphate of DMPC and the OmpT sidechains of the proposed LPS binding site, of the nearby Lys177 and of a second possible LPS site identified by examination of snapshots from the simulations. The results of this analysis are summarised in Table 3 and some snapshots are shown in Fig. 8B. From the table we can see that there are long-lasting (i.e. >8% of the duration of a given simulation) H-bonds from the DMPC molecules to all of the residues in the proposed LPS binding site in at least one simulation. The analysis suggests that the proposed site may also have an affinity for the phosphate part of the PC headgroup. The simulations also suggest that Lys177 may play a role in this binding site.

Interestingly, residues Tyr134, Glu136, Arg138, Arg175 and Lys177 form interactions with the detergent molecules (β-octylglucoside) present in the OmpT crystal structure. This would suggest that this region of the OmpT surface may have a general affinity for lipid and detergent headgroups. In the intracellular interfacial region the simulation results suggest that two aromatic sidechains (Tyr126 and Tyr189) may form a second lipid headgroup interaction site.

A more detailed analysis of the lipid/protein interactions in simulations OMPT1 and OMPT2 at the proposed LPS site shows little dynamics, with on average one DMPC molecule interacting via the phosphate moiety with each of the Arg138, Arg175 and Lys226 sidechains. It would be interesting to examine whether the diphosphate of Lipid A would form even more stable interactions, possibly with more than one basic sidechain at a time, thus supporting the proposed preferential binding of LPS at this site.

Given the importance of aromatic sidechains in membrane protein interactions with bilayers, we have also analysed the interactions of the two aromatic bands on the surface of OmpT with the phospholipid headgroups. The results for simulations OMPT1 and OMPT2 (Fig. 9) are broadly similar. In both cases, there are rather more aromatic/headgroup interactions at the periplasmic interface than at the extracellular. This analysis also reveals fluctuations in the number of such interactions on a ~2 ns timescale, reinforcing the conclusion above from analysis of numbers of H-bonds.

Overall, these simulations are encouraging in suggesting that 10 ns simulations can reveal key aspects of lipid-protein interactions. Of course, we must remember that the sites in the MD simulations are for interaction with a simple (phosphatidylcholine) bilayer, rather than the more complex lipid environment in the in vivo outer membrane. However, the simulation results do enable the formulation of experimentally testable hypotheses as to likely lipid interaction sites on the surface of OmpT.

Discussion

Biological Implications

The results of these simulations appear to support the novel mechanism for OmpT proposed by [START_REF] Kramer | Identification of essential acidic residues of outer membrane protease OmpT supports a novel active site[END_REF]. In particular, they demonstrate that even in the absence of bound peptide, the H212 and D83 sidechains bind a water molecule in an orientation that would favour nucleophilic attack on the scissile peptide bond. Simple peptide docking calculations indicate that a peptide can readily bind in the appropriate orientation at the active site for such a mechanism.

The simulations also reveal the role of water within the b-barrel. Although such water may well play a 'structural' role in stabilizing the transmembrane b-barrel fold, it also seems to be relatively mobile, that is the OmpT barrel forms a water-permeable pore. It is not clear whether or not this has a functional role. It is conceivable that the mobile water within the pore may play a role in 'delivering' water molecules to/from the active site at the external end of the barrel. Alternatively, it may be that the water is playing mainly a structural role and the water permeability of the OmpT barrel is simply tolerated (in an evolutionary sense) given the high permeability of the outer membrane due to the presence of porins. It would be possible to design suitable 'pore-perturbing' mutants that might enable these two alternatives to be tested experimentally.

In terms of interaction of OmpT with its membrane environment, the simulations are quite revelatory. In particular, although the simulations have been performed in a relatively simple DMPC bilayer, they appear to reveal the presence of specific lipid interaction sites on the surface of the OmpT barrel. These sites include the one proposed [START_REF] Vandeputte-Rutten | Crystal structure of the outer membrane protease OmpT from Escherichia coli suggests a novel catalytic site[END_REF], on the basis of structural homology with FhuA [START_REF] Ferguson | A conserved structural motif for lipopolysaccharide recognition by procaryotic and eucaryotic proteins[END_REF], to form a site for LPS binding. This result is of interest in terms of revealing the ability of extended MD simulations to provide meaningful information on protein/lipid interactions (Domene et al., 2003b).

Simulation Methodology and Future Directions

From a methodological standpoint these simulations reveal some small differences in behaviour according to whether long range electrostatic interactions are approximated via a simple cutoff or treated more accurately via the PME method. (A more detailed comparison of the effects of this aspect of the methodology on outer membrane protein simulations will be presented elsewhere). However, we note that although more accurate than the cutoff treatment, the PME method is not without possible artefacts [START_REF] Weber | Molecular dynamics simulations of a polyalanine octapeptide under Ewald boundary conditions: Influence of artificial periodicity on peptide conformation[END_REF], especially in the context of membrane protein simulations [START_REF] Bostick | The implementation of slab geometry for membranechannel molecular dynamics simulations[END_REF], which makes some more methodological development in this area desirable.

Our simulation results suggest, in terms of the tilting of the protein relative to the bilayer, that some degree of lipid/protein mismatch may occur for OmpT in DMPC bilayers. More extended simulations are required to obtain a more statistically significant picture of this possible mismatch, and also to explore how the local bilayer geometry may adjust to it. This information could be used to predict experimentally verifiable distances between lipid headgroups and specific protein residues. Given the relaxation time of the protein position in the bilayer this is currently beyond our computational means.

Both this mismatch, and the observation of specific protein/lipid interactions hastens the need for simulations of outer membrane proteins in a more realistic model of the outer membrane than a simple phosphatidylcholine bilayer. Initial modelling and simulation studies of LPS bilayers have been performed [START_REF] Katowsky | Molecular modelling of bacterial deep rough mutant lipopolysaccharide of Escherichia coli[END_REF][START_REF] Lins | Computer simulation of the rough lipopolysaccharide membrane of Pseudomonas aeruginosa[END_REF] and simulations of a more realistic OM model should be feasible in the near future.

The simulations reported here provide valuable clues as to the catalytic mechanism of OmpT. However, in order to characterise this more fully it will be necessary to perform simulations using QM/MM methods [START_REF] Mulholland | Ab initio QM/MM study of the citrate synthase mechanism. A low-barrier hydrogen bond is not involved[END_REF][START_REF] Ridder | Modeling biotransformation reactions by combined quantum mechanical/molecular mechanical approaches: from structure to activity[END_REF]. A first approximation to this may be to perform molecular mechanics based simulations that allow for dynamic protonation/deprotonation of water and sidechains [START_REF] Wu | Computer simulation studies of proton transport in a synthetic ion channel[END_REF]. Table 1: The values indicate percentages of total simulation time for which a given H-bond was present.

Tables

An interaction is recorded only if it was present for >30% in at least one simulation. C Dominant secondary structure analysis for simulations OMPT1, OMPT2 and OMPT3

showing residues that maintain their secondary structure as determined via DSSP [START_REF] Kabsch | Dictionary of protein secondary structure:-pattern-recognition of hydrogen-bonded and geometrical features[END_REF] for more than 90% of the total simulation time. Red = b-strand, green = bend, yellow = turn, light grey = coil. The b-strands and extracellular loops are labelled. Results of docking calculations, showing a selected OmpT-ARRA complex in which interacting protein residues are in light blue, the peptide is in green, and the scissile peptide bond is highlighted in yellow. D83 is protonated and a water molecule between D83 and H212 is shown.

The putative catalytic mechanism is indicated via arrows. The I170 residue is not shown for clarity, but its approximate sidechain position is indicated by an asterisk. The number of OmpT/lipid H-bonds vs. time for simulations OMPT1 (black line) and OMPT2

(grey line).

Figure 8:

A Comparison of the covalent structures of DMPC vs. lipid A (the latter after [START_REF] Ferguson | A conserved structural motif for lipopolysaccharide recognition by procaryotic and eucaryotic proteins[END_REF]).

B Snapshots showing the interactions between basic residues and DMPC at the putative lipid A binding site for the three residues Lys226, Arg175 and Arg138.

A

  OmpT fold showing: (i) the four active site residues (i.e. Asp83, Asp85 on the left, and Asp210, His212 on the right) in red; (ii) the basic lipid binding sites formed by Arg and Lys residues in orange; and (iii) the upper and lower amphipathic aromatic belts formed by Trp andTyr residues in green. The position of the lipid bilayer is delineated by horizontal lines, with the extracellular space located at the top and the periplasm at the bottom. B A schematic crosssection of OmpT indicating the approximate position of a potential tetrapeptide substrate (in black). C Snapshots from the OMPT2 simulation at t = 0 and 10 ns. The protein is blue (cartoon representation), the lipid bilayer polar headgroups are shown as red and orange spheres, the hydrophobic tails as green lines, water molecules as blue dots, and Na + and Cl -ions as purple and cyan spheres respectively. Three regions are indicated: w = water, i = interface; h = hydrophobic core.

Figure 2 :

 2 Figure 2:

Figure 3 :AB

 3 Figure 3:

Figure 4 :ABC

 4 Figure 4:

  Figure 5:

Figure 6 :AC

 6 Figure 6:

  Figure 7:

Figure 9 :

 9 Figure 9:

  Summary of Simulations

	Name	PME or Cutoff	Duration (ns)	Cα RMSD* (Å)
				b-strands	Not b-strands
	OMPT1	PME	10	1.00 ± 0.07	2.61 ± 0.24
	OMPT2	Cutoff	10	1.20 ± 0.08	3.58 ± 0.28
	OMPT3a	Cutoff	2.5	0.72 ± 0.05	1.85 ± 0.32
	OMPT3b	Cutoff	2.5	0.77 ± 0.06	2.12 ± 0.31
	OMPT3c	Cutoff	2.5	0.94 ± 0.15	2.05 ± 0.43
	*The Cα RMSD is calculated, relative to the starting structure, for the final 2 ns of each
	simulation. Each simulation system consisted of 63205 atoms, corresponding to monomeric
	OmpT, 19 water molecules taken from the X-ray structure, 33 Na + ions, 30 Cl -ions, 249 lipid
	(DMPC) molecules and 16191 waters in total.

Table 2 :

 2 H-bond Interactions at the Active Site

	Interaction	OMPT1	OMPT2	OMPT3a	OMPT3b	OMPT3c
	D210-H212	100	100	100	100	100
	D210-E27	-	-	1.6	40.5	-
	D83-S99	7.0	31.9	2.4	13.5	31.7
	D83-H101	-	31.5	2.4	1.6	0.8
	D83-w-H212	15.8	52.1	78.6	67.5	49.2
	D83-w-D85	72.9	58.1	93.7	81.7	96.0
	D85-S40	81.4	55.7	77.0	93.7	98.4

Table 3 :

 3 Dominant Lipid/Protein Interactions

	Protein-lipid interactions (given as percentage of total simulation time for which an H-bonding
	interaction is formed between the phosphate of DMPC and a protein sidechain) are listed for: (i)
	'proposed site' = residues identified by analogy with the Fhua LPS binding site; (ii) 'vicinity' =
	interactions with Lys177 which sits near to the proposed site; (iii) '2 nd site' = a putative second
	LPS-binding site; and (iv) 'lower belt' = interactions with the lower (periplasmic) half of the
	protein.

Baaden & Sansom, Fig. 8

Acknowledgements.

The authors thank José Faraldo-Gómez for providing the water bridge analysis program and Sundeep Deol for assistance in analysing lipid-protein interactions. MB thanks the EU for a grant (Contract-N o QLK2- CT-2000-51210). Work in MSPS's laboratory is supported by grants from the Wellcome Trust, the BBSRC and the EPSRC. We thank the Oxford Supercomputing Centre for access to resources.